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ABSTRACT 

A  least  squares  algorithm  is  developed  to  solve  for  the 
trajectory  and  transponder  array  coordinates  of  the  current 
velocity  profiler,  Pegasus.  Measurement  residuals  and 
parameter  precision  are  computed  for  data  quality  analysis. 
Travel  times  from  a  maximum  of  four  seafloor  transponders, 
pressure  sensor  depths,  and  transponder  positions  are  input 
with  their  respective  accuracy  estimates.  The  algorithm  is 
used  to  analyze  a  2250  m  profile  from  the  Monterey  Canyon  with 
four  transponders,  one  of  which  had  not  been  positioned.  This 
transponder's  unknown  position  is  found  and  problems  in  the 
other  array  coordinates  identified.  Transponder  coordinate 
precision  improves  by  factors  of  ten  in  the  horizontal  and 
five  in  depth,  to  about  13  m  (Drms)  and  2  m  (la)  respectively. 
Trajectory  precision  is  about  7  m  (Drms)  horizontal,  with  high 
correlation  between  points.  Thus,  the  precision  of  horizontal 
velocity  components,  determined  by  time  differencing  points, 
is  better  than  11  cm/s  (la) .  Depth  precision  is  better  than 
3m  (lo) ,  except  in  the  deeper  portions  where  anomalous 
pressure  residuals  near  the  depth  of  the  transponder  array 
suggest  systematic  pressure  errors  needing  further  study. 
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I.   INTRODUCTION 

The  Oceanography  Department  of  the  Naval  Postgraduate 
School  is  presently  conducting  a  study  of  ocean  current 
velocities  in  the  waters  off  the  central  California  coast  in 
the  area  of  the  Monterey  Canyon.  During  the  course  of  the 
study,  several  stations  have  been  established  for  the 
collection  of  data.  Station  10,  northernmost  in  this  chain  of 
stations  and  the  primary  focus  for  analyses  presented  in  this 
thesis,  contains  four  seafloor  transponders,  one  of  which  was 
considered  to  have  an  essentially  unknown  position  at  the 
beginning  of  the  work  described  here. 

The  main  oceanographic  instrument  used  in  this  survey  is 
the  dropsonde  Pegasus.  As  it  descends  and  ascends  through  the 
water  column,  it  records  at  16  second  intervals,  both  the  data 
from  oceanographic  sensors  and  the  time  for  return  travel  of 
acoustic  signals  from  transponders  previously  set  in  place  on 
the  ocean  floor. 

In  the  past,  Pegasus-generated  velocity  studies  have  used 
the  raw  data  in  various  ways.  The  trajectory  of  the 
instrument  has  been  fitted  to  a  curve  which  was  then  differen- 
tiated with  respect  to  time  in  order  to  produce  the  horizontal 
velocity  components  (Halkin  et  al.,  1985).  Alternatively,  and 
as   is  the  present  practice  at  NPS,  velocities  have  been 


calculated  from  position  differences  and  then  smoothed  with  a 
spatial  filter  such  as  a  seven  point  running  average. 

This  study  outlines  a  method  of  determining  the  current 
velocities  and  their  precision  by  applying  a  least  squares 
adjustment  procedure,  adapted  from  geodetic  survey  techniques, 
whereby  all  observational  data  are  used  simultaneously. 
Specifically,  three  basic  problems  will  be  analyzed: 

*  The  precision  of  the  transponder  coordinates  from  which 
the  position  of  Pegasus  is  derived. 

*  The  precision  of  the  Pegasus  positions  and  consequently 
the  current  velocities  derived  from  Pegasus  navigation. 

*  The  depth  at  which  pressure  becomes  critical  in  solving 
for  velocities. 

Discussions  of  this  topic  will  appear  in  the  following 

sequence: 

*  Chapter  II  provides  background  information  on  the 
transponder  network  and  the  Pegasus  system  and  its 
positioning. 

*  Chapter  III  discusses  sources  of  error  resulting  from 
systematic  and  observational  inaccuracies  in  the  data. 

*  Chapter  IV  gives  a  brief  explanation  of  the  least  squares 
adjustment  process. 

*  Chapter  V  describes  the  final  procedure  used. 

*  Chapter  VI  discusses  the  results  obtained. 

*  Chapter  VII  provides  conclusions  and  recommendations  for 
future  consideration. 

*  Appendices  contain  algorithms  and  explanations  of  the 
Fortran  programs  used. 


II.   BACKGROUND 

The  area  being  surveyed  is  composed  of  ten  stations 
located  just  off  the  California  coast  approximately  between 
36° 05'  and  36° 39'  North  Latitude,  and  122° 08'  and  124° 13'  West 
Longitude.  Of  particular  interest  here  is  Station  CIO, 
northernmost  in  the  chain  and  located  in  the  Monterey  Canyon 
about  ten  miles  west  of  Monterey  (see  Figure  1) . 

A.   TRANSPONDER  NETWORK 

At  this  station  four  transponders  were  deployed,  horizon- 
tally separated  by  distances  somewhat  equivalent  to  transpon- 
der depths  (around  2000  m)  .  Three  of  these  transponders, 
responding  at  11.5  kHz,  12.0  kHz,  and  12.5  kHz,  respectively, 
were  located  by  the  survey  vessel  according  to  traditional 
methods  as  described  below.  However,  the  bandwidth  of  the 
shipboard  receiving  instrument  was  too  narrow  to  pick  up  the 
13.5  kHz  signal  from  the  fourth  transponder,  and  therefore  its 
position  was  largely  unknown.  (Table  1  defines  transponder 
abbreviations,  their  frequencies,  and  their  approximate 
depths.   Figure  2  shows  the  survey  net  for  Station  C10.) 
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Figure  1.   General  Area  of  Station  CIO 


TABLE    1 
CIO    TRANSPONDERS 


Transponder 

1 

Tl 

12.5 

kHz 

1880 

m 

Transponder 
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T2 

12  .  0 

kHz 

1990 

m 

Transponder 

3 

T3 

11.5 

kHz 

2230 

m 

Transponder 

4 

T4 

13.5 

kHz 

1890 

m 

*  TRANSPONDERS 
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Figure  2.   Station  CIO  Survey  Net 


B.   SURVEYING  A  TRANSPONDER  NETWORK 
1 .   Transponder  Depth 

To  determine  each  available  transponder  depth,  the 
survey  ship  homes  in  on  the  instrument's  signal  until  a 
minimum  time  is  recorded,  thereby  assuming  the  ship  to  be 
directly  over  the  instrument  (see  Figure  3).  An  appropriate 
harmonic  mean  sound  velocity  multiplied  by  one  way  signal 
travel  time  produces  the  depth  of  the  transponder. 
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Figure    3.       Determining   Transponder   Depth 


2 .       Baselines 

To  establish  the  relative  horizontal  x  and  y  coordi- 
nates of  transponders,  baselines  are  determined  between  the 
instruments.  Running  at  a  relatively  slow  speed  (say  4-6 
kts) ,  the  survey  vessel  attempts  to  cross  each  baseline  at  a 


90°  angle.  The  process  is  repeated  in  the  opposite  direction. 
This  procedure  is  then  duplicated,  several  more  times  if 
possible.  Careful  observation  of  the  analog  trace  that 
records  round-trip  travel  times  of  the  acoustic  signal  between 
the  survey  vessel  and  the  two  transponders  indicates  when  the 
shortest  distance  has  been  reached,  as  shown  in  Figure  4. 
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Figure  4.   Analog  Trace  of  a  Baseline  Crossing 


Signals  show  up  on  the  trace  as  parabolas  approaching 
one  another.  They  will  lie  vertically  in  line  if  the  baseline 
is  crossed  at  a  90°  angle  and  therefore  the  shortest  distance  to 
each  of  them  is  reached  simultaneously.  At  this  point  the 
ship  will  lie  in  a  vertical  plane  containing  the  two 
transponders.   The  sum  of  the  two  horizontal  ranges  will   be 


at  a  minimum  and  therefore  establish  the  baseline  (see  Figure 
5)  . 


Figure  5.   Baseline  Determination 

The  time  of  crossing,  ship's  course,  latitude/ 
longitude,  and  Loran  coordinates  are  noted.  Time  units  for 
each  minimum  range  are  measured  and  then  multiplied  by  the 
appropriate  harmonic  mean  sound  velocity  to  obtain  slope 
distances.  These  distances,  the  depths  of  the  two  transpon- 
ders, and  the  ship's  course  are  then  used  to  calculate  the 
horizontal  baseline  length  and  the  azimuth  between  the  two 
transponders.  Finally,  a  local  xy  plane  coordinate  system  is 
established,  setting  the  position  of  one  of  the  transponders 
at  0,0,  and  positioning  the  other  transponder  relative  to  it. 

Many  factors  combine  to  make  the  information  on  this 
analog  trace  subjective  and  inaccurate.   The  analog  paper  may 
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be  set  to  progress  at  different  rates  producing  sweeps  of 
0.5  s,  1  s,  or  2  s,  etc.  The  trace  may  actually  "wrap  around" 
one  or  more  times  as  the  time  interval  gets  longer  and  longer, 
making  page  integer  resolution  difficult.  At  a  one  second 
sweep,  the  full  page  represents  (after  conversion  from  time  to 
distance)  meters  traveled  in  one  second,  i.e.,  1500  m/s. 

Furthermore,  if  the  baseline  is  not  crossed  perpendic- 
ularly, a  time  offset  between  minimum  ranges  is  observed  which 
decreases  the  length  of  the  baseline.  This  error  can  be 
removed  to  some  extent  by  multiplying  the  time  offset  by  the 
ship's  speed  and  then  by  geometric  relationships,  computing 
the  correction.  This  situation  is  illustrated  in  Figure  6 
where  T1-T2  represents  the  actual  baseline,  A,  B  represent 
respective  points  on  the  analog  trace  where  minimum  slant 
ranges  are  indicated,  and  a,  b  represent  respective  minimum 
slant  ranges  to  each  transponder. 
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Figure  6.   Time  Offset 


In  addition  to  the  above  possibilities  for  error,  an 
inaccurate  course  heading  at  the  time  of  baseline  crossing 
will  alter  the  geometry  of  the  array,  resulting  in  inaccurate 
coordinates.  Further  discussions  of  transponder  depth  and 
baseline  determinations  can  be  found  in  Kuo  (1985) ,  McKeown 
(1975) ,  and  Hart  (1967) .  For  an  assessment  of  transponder 
network  accuracy,  refer  to  Chapter  III. 

Figure  7  illustrates  the  CIO  transponder  network  and 
the  approximate  position  of  Pegasus  drop  #157  positioned  in  a 
local  coordinate  system  using  Tl  (12.5  kHz)  as  0,0. 
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Figure  7.   The  CIO  Local  Coordinate  System 
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C.   DESCRIPTION  OF  THE  PEGASUS  SYSTEM 

The  oceanographic  instrument  Pegasus  is  a  free-falling, 
acoustically-tracked  velocity  profiler.  It  was  developed  by 
H.T.  Rossby  and  D.  Dorson  at  the  University  of  Rhode  Island  a 
decade  ago  to  fill  a  need  for  an  instrument  which  would  accu- 
rately record  the  fine  scale  vertical  structure  in  the  ocean 
and  prove  inexpensive  and  easy  to  handle  at  sea.  It  provides 
a  means  to  measure  absolute  velocity  components  throughout  the 
water  column.  (For  the  seminal  paper  on  Pegasus,  see  Spain  et 
al.,  1981.)  Prior  uses  of  Pegasus  in  studies  of  the  Gulf 
Stream  off  the  east  coast  of  the  U.S.  have  been  documented  in 
a  paper  and  data  report  by  Halkin  et  al.  (1985). 

The  position  of  the  instrument  as  it  descends  and  ascends 
is  determined  by  travel  time  of  signals  emitted  by  Pegasus  at 
10  kHz  every  16  seconds  which  are  heard  and  then  answered  at 
other  freguencies  by  transponders  previously  set  in  place  on 
the  ocean  floor,  and  by  pressure  readings.  These  time 
intervals,  and  temperature,  conductivity,  and  pressure  data 
from  oceanographic  sensors,  are  recorded  and  stored  in  a 
microprocessor-controlled  memory  in  the  instrument  to  be  later 
down-loaded  aboard  ship  when  the  Pegasus  is  recovered. 

The  Pegasus  model  currently  adapted  for  use  by  the  Naval 
Postgraduate  School  is  manufactured  by  Benthos.  It  is  housed 
in  a  17-inch  glass  sphere  protected  by  a  hard  hat  and, 
according  to  manufacturer  claims  (Benthos,  1989) ,  provides 
operation   to   full   ocean   depth,   integral   flotation,   no 
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corrosion,  and  large  battery  capacity  allowing  for  100  deploy- 
ments without  opening.  (A  normal  NPS  school  cruise  uses  about 
20  deployments.)  It  weighs  approximately  45  kg  in  air  and 
carries  eight  kg  of  expendable  weight  which  sink  the  profiler 
in  this  present  study  at  about  27  m/min.  Six  receiver 
channels  for  communication  with  transponders  are  available 
with  frequencies  located  at  0.5  kHz  intervals  ranging  from 
11.5  kHz  to  14.0  kHz.  NPS  supplied  SEA  BIRD  model  SBE-3 
temperature  sensor  and  SEA  BIRD  model  SBE-4  conductivity 
sensor  are  externally  mounted  with  electrical  interfaces, 
along  with  a  Paroscientif ic  model  410KT  pressure 
transducer  with  a  companion  Intelligent  Transmitter  circuit 
board. 

D.   PEGASUS  POSITIONING 

1 .   Computation  of  Pegasus  Coordinates 

The  position  of  Pegasus  at  each  16  second  interval  can 
be  determined  by  solving  a  system  of  equations  using  the 
formula: 

Time  =  Distance/Sound  Velocity 

Using  one-way  signal  times  from  two  transponders  (see 
Figure  8)  ,  the  depth  of  Pegasus  as  derived  from  pressure 
information,  and  an  appropriate  sound  velocity,  there  will 
then  be  two  equations  in  two  unknowns,  allowing  for  a  solution 
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PEGASUS 


Figure  8.   Determining  Pegasus  Position 

of  XP  and  YP  at  a  given  Pegasus  position.   The  solution  is 
derived  from  the  equations: 


Time!  =  [  (XP  -  Xx)£   +  (YP  -  Yx)z  +  (ZP  -  Zi)  ]  ' /V 


Time2   =    [  (XP    -    X2)2   +    (YP    -    Y2)2   +    (ZP    -    Z2)2]1/2/V 


where 


Timei    =  one-way  travel  time  from  transponders  1,2; 

V       =  effective  sound  velocity  along  the  path  between 
Pegasus  and  the  transponder; 

XP,YP,ZP  =  Pegasus  position  coordinates  (ZP  is  assumed 
known  from  the  simultaneous  pressure 
observations) ; 
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X1,Y1,Zi    =  transponder  coordinates  (assumed  known). 

The  traditional  oceanographic  method  outlined  above 
calculates  Pegasus  position  coordinates.  However,  it  lacks  a 
means  of  obtaining  precision  estimates  for  the  derived 
positions,  and  it  is  not  able  to  use  all  the  observed  data 
simultaneously.  Furthermore,  there  is  no  redundancy  in  the 
computational  process  and  thus  systematic  errors  or  blunders 
can  escape  unnoticed. 

2 .   Estimating  Horizontal  Current  Velocity 

Current  velocity  vectors  u  and  v  (describing  velocity 
in  the  x  and  y  directions,  respectively)  between  two  Pegasus 
positions,  can  be  obtained  by  taking  the  difference  between 
coordinates  and  dividing  by  the  time  interval: 

u  =  (XP2  -  XPJ/dt 

v  =  (YP2  -  YPJ/dt 

dt  =  16  s 

This  method,  like  the  one  mentioned  in  the  previous 
section,  does  not  provide  estimates  of  precision  and  does  not 
make  use  of  all  available  data. 
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III.   SOURCES  OF  ERROR 

Davis  et  al.  (1981,  pp.  15-20)  defines  measurement  error 
as  the  difference  between  a  measured  and  a  "true11  value  and 
describes  these  errors  as  being  generally  of  three  types: 

*  blunders  or  mistakes. 

*  systematic  errors. 

*  random  errors. 

Blunders  can  be  caused  by  carelessness,  equipment  failure, 
or  false  interpretation.  They  are  usually  large  enough  to  be 
easily  spotted  when  results  are  analyzed. 

Systematic  errors  follow  a  defined  pattern  or  system, 
which  when  discovered,  can  usually  be  expressed 
mathematically.  Such  factors  as  observer  limitations, 
instrument  imperfections  or  inadequate  calibration, 
meteorological  conditions,  or  poor  choice  of  mathematical 
model  can  produce  the  pattern  which  will  remain  consistent  as 
long  as  the  elements  of  the  system  remain  the  same. 
Systematic  errors  are  not  eliminated  by  repetition  of 
measurements.  They  must  be  ferreted  out  and  either  removed 
from  the  observations  or  their  effects  added  to  the  model. 
Thus,  to  begin  such  an  analysis,  serious  consideration  has  to 
be  given  to  locating  and  defining  errors  resulting  from 
systematic  and  observational  inaccuracies  in  the  data. 
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Random  errors  are  those  that  remain  after  mistakes  and 
systematic  errors  have  been  accounted  for.  The  values  of 
these  errors  should  be  unbiased  and  will  ideally  distribute 
themselves  about  a  mean  of  zero.  It  is  these  random  errors 
that  the  least  squares  process  attempts  to  minimize.  In  the 
analysis  described  in  this  thesis,  a  considerable  amount  of 
time  was  spent  in  determining  what  the  estimates  of  the  a  priori 

standard  deviations  of  these  random  errors  should  be,  both 
from  the  point  of  view  of  reasonable  physical  reality  and 
meaningful  results.  The  adjustment  technique  uses  these  a  priori 

estimates  to  weight  the  contribution  of  each  observation  to 
the  final  solution. 

The  following  potential  error  sources  were  studied  and, 
where  applicable,  appropriate  standard  deviations  of  the 
errors  were  then  entered  as  weights  into  the  least  squares 
adjustment  model: 

*  Transponder  coordinates. 

*  Signal  time. 

*  Velocity  of  sound. 

*  Pressure/Depth  relationship. 

A.   TRANSPONDER  SURVEY 

Of  prime  importance  is  the  precision  of  the  transponder 
positions.  Determining  these  coordinates  is  a  difficult  and 
time-consuming  process.  As  noted  before,  in  the  CIO 
transponder  net,  the  location  of  the  4th  transponder  could  not 
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be  determined  by  traditional  methods  because  the  band  width  of 
the  ship's  receiver  was  not  wide  enough  to  pick  up  the 
13.5  kHz  signal.  Also,  it  was  not  possible  at  the  time  to 
determine  more  than  two  of  the  three  baselines  between  the 
three  "known"  points,  thereby  precluding  a  mathematical 
closure  of  the  triangle. 
1 .   Transponder  Depth 

To  establish  a  reasonable  estimate  for  the  standard 
deviation  of  transponder  depth  error,  it  was  noted  that  any 
horizontal  offset  resulting  from  the  ship's  not  being  directly 
overhead  always  produces  a  positive  depth  error,  i.e.,  the 
slant  range  must  be  longer  than  the  vertical.  Figure  9  is 
similar  to  Spain  et  al.  (1981,  p.  1557)  and  illustrates  this 
process. 

The  solution  is  as  follows: 


(H  +  h)z  =  xz  +  H 


2Hh  +  h2  =  x2 


2         2 
_  =  h  +  _=  h(1  +  __, 


where: 


H  =  true  transponder  depth; 
h  =  depth  error. 


17 


y^ 


v^ 


h  = 


2H 


Figure  9.   Transponder  Depth  Error 


If  h  <<  H,  the  last  term  can  be  disregarded,  and 


n  ~  — 
~  2H 


A  realistic  estimate  of  the  ability  of  a  survey  vessel 
to  cruise  directly  over  a  transponder  by  using  this  method  is 
about  100  meters  (Schnebele,  1990a) .  Given  the  uneven  canyon 
terrain  of  this  survey  which  makes  resolving  the  position  of 
the  shoalest  depth  even  more  difficult,  it  seemed  reasonable 
to  ascribe  an  error  of  200  m  to  the  offset.  Thus,  if  the  true 
depth  lies  at  2000  m  and  there  is  a  horizontal  offset  of 
200  m,  there  will  be  a  depth  error  of  about  10  m. 
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Therefore,  a  oz  =  10  meters  was  chosen  for  transpon- 
ders 1,  2,  and  3.  Since  there  was  no  observed  depth  for 
Transponder  4,  a  depth  of  1900  meters  was  read  from  the  latest 
NOAA  chart  using  the  position  of  transponder  drop1  and,  since 
it  had  not  been  surveyed  in,  a  larger  standard  deviation  of 
100  meters  attached  to  it.2 


°(Tiz,T2z,T3z)  —  10  m 


oT<,z  =  100  m 


2 .   Baselines 

In  determining  baselines,  three  sources  of  error 
predominate: 

*  Minimum  range  reading  taken  when  the  ship  is  offset  from 
the  actual  baseline. 

*  Transponder  depth  error  propagating  into  the  computed 
baseline  length. 

*  Error  in  ship's  course  affecting  baseline  azimuth. 

In  a  simplified  ideal  situation  where  the  baseline  is 
crossed  at  a  90°  angle  in  the  center  of  the  line  and  the 
depths  of  the  two  transponders  are  the  same,  the  expected 


xThe  transponder  position  was  calculated  from  Loran 
signals  which  are  based  on  North  American  Datum  27.  The  NOAA 
chart  uses  NAD83.  Positions  may  vary  up  to  100  meters  between 
the  two  datums  in  this  locality. 

2After  the  October  17,  1989  earthguake,  a  depth  measure- 
ment of  1931  meters  was  made  of  the  4th  transponder  by  NPS 
scientist  Tarry  Rago  in  a  submersible. 
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error  in  baseline  length  from  measured  slant  ranges  can  be 
calculated  as  follows  (see  Figures  10  and  11)  . 
a.   Baseline  Error  Due  to  R 


Figure  10.   Baseline  Error 


B 
b 


=  distance  off  baseline  when  slant  ranges  were 
measured; 

=  true  baseline  between  transponders  (B  =  B1   +  B2) ; 
=  baseline  error  (b  =  b1  +   b2) ; 
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2     _     t,2 


(Bj    +    bj)'    =    Rz    +    B2 


2     _     T52 


28^!    +    b/    =    R 


-     d2 


b1(2B1    +    bj     =    R 


b  R2  R2  R2  bl  bl2 

Di  "  2b,  +b. e:     >2B:tl"    217    +(2b7)    '  •  •  • J 

2Bl(l    +  ^) 


~bl 
Since    (j^-  +    .  .  .)    <<    1, 

bx    ^ and,    similarly,    b2     ^ 


-  2B1  '  J '      *     -   2B2 


If,    as    stated   above,    Bx   =    B2   =    B/2,    then   bx   =   b2   =   b/2 


2  2  2  2 

.        ,.  JT  JT  R_,_L_         _L%  R_/ix 

Dl      D2    -   2B  2B2    -      2  ^B/2         B/2J    -     2  V 

b      2r2 

bR  =  — 


This  error  increases  the  length  of  the  baseline. 
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b.   Baseline  Error  Due  to  h 


* 
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Figure  11.   Relationship  of  Baseline  and  Depth  Errors 

H   =   true  transponder  depth; 

h    =   depth  error; 

Sa   =   slant  range  from  ship  to  transponders  1,2. 


si2  =  Hi  +  Bi2 


(1) 


Sr  =  (Hx  +  h^  +  (Bi  +  bx) 


(2) 


Sx  remains  constant. 

Eg.  (2)  -  Eg.  (1) : 
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But 


0    =    2H1h1    +    hx2    +    2B!bi    +    b/ 


2Blbl(1   +2BL)     =    ^VlC1   +^±) 
1  2H 


Blbl 


H   h 


(1+2H7} 


(1    +2B7) 


B1b1  h1  h1    2 

p— -  -    B1b1(l    -    j—   +     <2h~)       ~    •••) 

(1   + — — ) 
v         21^' 

-H    h  b  b       2 

E7 Hlhl(1  "  2B7+    <2b7>     "   •••> 

(1  +2B7' 
b  L  h 

Since    (__L+    ...)    «    1,    and    ( —i.  +    ...)    «    l,    then 

2B  2H, 


B^!   *,   -Hihi 


bi     * 


Hlhl 
B, 


If    Hi    =    H2    =    H,    h,    =    h2    =    h,     Bj    =    B2    =    B/2,    ba    =    b2    =    b/2,    th 


en 


bh   =   bx   +   b2 


Bl         B2 


-Hh(=r~r  + 


B/2         B/2 


bh   ~    - 


4Hh 
B 
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As  shown  in  Figure  11,  the  hypotenuse  representing 
the  given  slant  range  of  the  signal  must  be  the  same  in  both 
triangles  ADE  and  ACF.  FC  must  equal  ED.  Therefore,  if  the 
assumed  depth  is  increased  by  error,  then  the  computed 
baseline  must  consequently  decrease, 
c.   Total  Baseline  Error 

Thus  the  total  baseline  error  bR  +  bh  becomes: 


blotal  =  bR  +  bh    ~  !(2R2   "   4Hh) 


Using  typical  figures  of:  B  =  2000  m,  H  =  2000  m, 
h  =  10  m,  and  R  =  100  m,  it  should  be  possible  to  obtain  a 
value  for  bTotal  =  ±30  m. 

If  the  transponder  depths  are  not  the  same,  and 
the  baseline  is  not  crossed  right  at  the  center  but  is  crossed 
on  a  perpendicular  heading,  the  baseline  error  will  be: 

2  Hh,         H0h_ 

.  ?_r_L       JL\         <    1   1    +      2    2) 

^Totai     ~      2     B1  B2;  K    B±  B2    ' 

In  addition  to  baseline  errors  discussed  above,  an 
incorrect  azimuth  will  add  even  more  uncertainty.  An  azimuth 
error  of  0.5°  over  a  baseline  of  2.0  km,  for  example,  will 
produce  an  error  in  position  of  one  end  of  the  baseline  with 
respect  to  the  other  of  18  meters. 
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Based  on  the  above  analysis,  it  was  initially 
assumed  that  the  relative  horizontal  coordinates  of  Tl,  T2 , 
and  T3  would  be  known  to  approximately  ±30  m.  Hart  (1967), 
for  example,  claims  that  with  careful  repetitive  procedures 
involving  all  baselines  in  a  transponder  array,  relative 
accuracies  of  the  order  of  five  meters  can  be  achieved. 

As  preliminary  data  analysis  proceeded  with  the 
CIO  network,  it  became  obvious  that  there  were  major  inconsis- 
tencies with  the  given  baseline  data.  Preliminary  adjustments 
(see  Chapter  V)  would  not  converge  when  the  given  transponder 
coordinates  were  used.  Further  analysis  suggested  both  length 
and  azimuth  problems  with  the  given  T2-T3  baseline.  The  fact 
that  the  T3-T1  baseline  had  never  been  determined  led  to  a 
situation  in  which  no  positional  redundancy  existed  in  the  Tl, 
T2 ,  and  T3  network.  As  a  result  of  these  problems,  initial 
positions  of  Tl,  T3 ,  and  T4  were  determined  from  the  Loran 
coordinates  of  their  drop  sites,  while  T2  was  computed 
relative  to  Tl  from  baseline  survey  data. 

As  a  consequence  of  this  less  than  ideal 
procedure,  it  was  considered  that  standard  deviations  of  ±100 
m  should  be  allocated  to  the  transponder  coordinates  of  T2 , 
T3 ,  and  T4 ,  while  Tl  would  be  assumed  fixed  in  order  to 
provide  a  positional  datum  (albeit  a  rather  arbitrary  one)  for 
the  array.  The  100  m  allocated  for  positional  standard 
deviations  seemed  reasonable  in  the  light  of  both  the  known 
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positional  accuracy  of  Loran  and  the  added  uncertainty  from 
drift  as  the  transponders  settled  to  the  bottom. 

In  summary,   therefore,   the  following  standard 
deviations  were  assumed: 


a  (T2    ,T3    ,T4    )  =  100  m 
x,y'   x,y'   x,y 


oT1     =  0.01  m 
x,y 


B.   TIMING 

The  time  observations  required  a  consideration  of  how  well 
the  acoustic  time  intervals  could  be  established  between 
Pegasus  and  the  transponders.  The  travel  times  were  assumed 
to  be  independent  of  one  another.  The  Pegasus  manual  states 
that  the  transponders  have  a  12.5  ms  output  pulse  delay 
accurate  to  within  ±0.5  ms  (Oceanographic  Instrument  Systems) . 
For  the  initial  transponder  adjustment,  oTime  =  0.002  s  was 
chosen  as  a  reasonable  weight  (see  Chapter  V) ,  and  then  later 
reduced  to  0.0005  s  for  the  full  run  with  improved  transponder 
coordinates. 


Olimel  =  0.002  s 


°Time2     =     0.0005     S 
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C.   VELOCITY  OF  SOUND 

Fofonoff  (1963)  states  that  velocity  of  sound  is  a 
function  of  the  thermodynamical  and  chemical  state  of  the 
water  and  is  determined  by  using  any  complete  set  of  variables 
of  state  such  as  temperature,  pressure,  and  salinity.  He  goes 
on  to  observe  that  the  precision  of  sound  velocity  must 
consider  not  only  that  of  the  oceanographic  data,  but  also  of 
the  empirical  formula  used  to  convert  these  measurements  to 
sound  speed. 

Pegasus  provides  temperature  and  conductivity  data  which 
is  converted  to  a  sound  velocity  profile  for  the  full  range  of 
the  drop  by  algorithms  published  in  the  Unesco  Technical 
Papers  in  Marine  Science  #44  (Fofonoff  and  Millard,  1983,  pp. 
11-12,  49) . 

In  dealing  with  slant  ranges,  as  in  this  case,  acoustic 
refraction  should  be  considered  as  well.  Spain  et  al.  (1981, 
pp.  1562-1564)  presents  an  equation  derived  by  Vaas  (1964) 
which  corrects  for  the  effects  of  ray  bending,  taking  into 
consideration  the  relative  positions  of  the  projector  and  the 
receiver  even  when  they  are  at  similar  depths.  This  equation 
is  stated  to  have  an  accuracy  of  better  than  0.2  m/s  for  all 
angles  and  depths. 

Another  excellent  source  for  acoustic  refraction  informa- 
tion in  a  survey  situation  is  the  SASS  Accuracy  Study 
Simplified  Ray  Bending  Correction  (General  Instrument  Corp. , 
1975)   and  its  follow-up  Ray  Bending  Correction  for  Depth 
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Sounders,  An  Informal  Approach  (General  Instrument  Corp., 
1976) .  These  two  papers  develop  the  equations  for  an 
effective  sound  velocity  in  detail  and  include  graphs  showing 
the  errors  relative  to  lateral  angles. 

In  this  study  it  was  found  that  the  use  of  an  average 
harmonic  mean  sound  velocity  produced  acceptable  results.  To 
obtain  this  average,  the  sound  velocity  profile  was  used  to 
calculate  harmonic  sound  velocity  profiles  with  respect  to 
each  transponder.  These  four  harmonic  mean  profiles  were  then 
averaged  (see  Figure  12)  and  the  resulting  profile  used  to 
calculate  the  approximate  X  and  Y  coordinate  of  each  Pegasus 
position.  The  harmonic  mean  profiles  differ  for  each 
transponder  because  of  the  wide  difference  in  deployment 
depths.  Harmonic  means  at  the  depths  found  here  differed  by 
about  1.0  m/s  from  the  average  used,  representing  a  bias  of 
less  than  0.7  m  in  a  typical  one  km  path  length. 

Refraction  was  not  regarded  as  significant  for  this  study 
because  of  the  relatively  short  path  lengths  and  small  sound 
velocity  gradients  encountered  (Schnebele,  1990b) .  In  deeper 
waters  with  longer  paths,  refraction  may  be  significant. 

The  harmonic  mean  used  in  this  model,  therefore,  assumes 
no  refraction.  The  small  error  that  it  introduces  was  felt  to 
be  inconsequential  compared  to  the  uncertainties  in 
transponder  coordinates  and  signal  travel  time  measurements. 
Thus  the  harmonic  means  used  in  the  adjustment  were  assumed  to 
be  exact  quantities  without  significant  random  error. 
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Figure  12.   Harmonic  Sound  Velocity 

D.   PRESSURE/DEPTH  RELATIONSHIP 

The  algorithm  used  to  convert  pressure  to  depth  is  the 
Saunders  and  Fofonoff  formula  which  uses  the  hydrostatic 
equation  and  the  Knudsen-Ekman  equation  of  state  (Fofonoff  and 
Millard,  1983,  p.  25).  The  formula  includes  variation  of 
gravity  with  latitude  and  depth  and  assumes  standard  sea 
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water.   This  formula  is  said  to  deviate  by  only  0.08  meters  at 

5000  decibars  from  estimates  based  on  EOS80,  a  considerably 

smaller  error  than  those  in  pressure  measurements  as  shown 

below. 

The  manufacturer  of  the  Pegasus  pressure  transducer  claims 

a  typical  accuracy  under  difficult  environmental  conditions  of 

0.02%  (Paroscientif ic,  1986).   At  our  full  scale  of  2200  m, 

this   would   compute   to   about   0.4  m.     As   the   analysis 

progressed,  further  inconsistencies  in  the  data  began  to  lead 

to  a  suspected  bias  in  pressure  measurements,  perhaps  caused 

by  a  temperature  hysteresis  as  Pegasus  drops  and  rises  through 

the  deepest  part  of  the  water  column.  Therefore,  a  o   =2.00 

pz 
m  was  taken  to  be  the  standard  deviation  of  the  Pegasus  depth 

measurements,  i.e., 


o    =2.00. 

pz 
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IV.   LEAST  SQUARES  ADJUSTMENT 

The  specific  task  of  this  study  was  essentially  to  analyze 
three  basic  problems: 

*  The  precision  of  the  transponder  coordinates,  especially 
those  of  T4 ,  from  which  the  position  of  Pegasus  is 
derived. 

*  The  precision  of  the  Pegasus  positions  and  consequently 
the  current  velocities  derived  from  Pegasus  navigation. 

*  The  depth  at  which  pressure  becomes  critical  in  solving 
for  velocities. 

To  provide  a  means  of  answering  these  questions,  a  Least 
Squares  Adjustment  Program  was  developed  by  Dr.  John  Hannah.3 
This  adjustment  program,  as  well  as  a  brief  explanation,  are 
provided  in  Appendix  A. 

Least  squares  adjustment  techniques  are  used  to  determine 
unique  solutions  for  unknown  parameters  when  there  are 
redundant  observations,  i.e.,  more  than  necessary  to  specify 
the  model  (Davis  et  al.,  1981,  pp.  38-39).  The  least  squares 
estimator  is  an  unbiased  minimum  variance  estimator  which  is 
unique,  mathematically  easy  to  derive,  and,  when  compared  to 
other  estimation  techniques,  leads  to  a  smaller  dispersion  of 
random  errors.   It  also  is  distribution  free  in  the  sense 


3Adjunct  Research  Professor,  CNOC  Chair  in  Mapping, 
Charting,  and  Geodesy,  U.S.  Naval  Postgraduate  School, 
Monterey,  CA. ,  1988-1990. 
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that  a  distribution  is  needed  only  for  confidence  interval 
testing  (Hannah,  1989) . 

The  mathematical  model  for  the  Pegasus  data  adjustment  has 
the  general  form: 

Lai  =  Fx(xa) 

La2  =  F2(xa) 

in  which  the  first  set  of  observations  comes  from  Pegasus 
itself  in  the  form  of  return  signal  travel  times  from  the 
transponders  to  the  instrument.  The  second  set  comes  from  an 
a  priori   knowledge  of  any  of  the  system  parameters  such  as  the 

positions  of  the  transponder  coordinates,  and  each  Pegasus 
depth  as  calculated  from  pressure.  In  general  terms,  any  set 
of  adjusted  observations  expressed  as  a  function  of  a  set  of 
adjusted  parameters  can  be  written  in  a  matrix  expression: 

La  =  F(xa) 

The  following  is  adapted  from  a  least  squares  adjustment 
procedure  written  by  Hannah  (1981)  and  is  used  here  with  the 
author's  permission. 


...These  functions  may  be  linearized  by  taking  a  first  order 
Taylor  series  expansion  about  some  approximate  values  for 
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the  parameters,  X.   The  first  function  then  becomes 

DF1 

**   +  vi  =  ~^r  i        (Xa-x0)  +  F(xo) 

a  X  =Xn 
a  0 


or 


A 


or 


A 

Vj  =  A2X  +  Lt 


in  which 


8F3         a 
Ai  =  — i  |      ,   X  =  (Xa-X0) 

a  X  =Xn 
a  0 

and 

Lx  =  L0  ~  Lb 

Similarly,  the  second  function  may  be  linearized  to  give 

a 
V2  =  A2X  +  L2 


In  the  above,  Vj  and  V2  are  the  vectors  of  residuals  arising 
from  the  observations  L^1  and  L^2  with  the  adjusted  values  of 
the  parameters  being  given  by  the  vector  Xa.  Assuming  that 
the  two  sets  of  observations  are  uncorrelated,  then  the 
least  squares  minimum  variance  estimate  of  X  based  on  these 
two  sets  of  observations  is  given  by  minimizing  the  function 

A  A 

(J)  =  V^P^  +  V2TP2V2  -  2K1T(V1-A1X-L1)  -  2K2T  (V2-A2X-L2) 

A 

with  respect  to  the  unknowns  V1#  V2,  Klr  K2,  and  X.  The 
weight  matrices  for  each  set  of  observations  are  given  by  Px 
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and  P2,  in  which  Px  =  E^-l  and  P2  =  E2  -1,  i.e.  the  inverse  of 

L  b  L  b 

the  variance-covariance  matrices  of  the  observation  sets. 
Infinitely  large  variances  are  applied  to  non-weighted 
parameters  resulting  in  zeros  in  the  corresponding  diagonal 
elements  of  P2. 

After  minimizing  <J>  and  eliminating  the  unknowns  Klf  K2/ 
Vj,  and  V2/  the  least  squares  estimate  for  X  is  given  by  the 
solution  of  the  normal  equations 

(A1TP1A1    +    A2TP2A2)X    +     (A^PiLi    +    A2TP2L2)     =    0 


Since,  however,  the  A2  matrix  arises  from  direct  observa- 
tions on  the  parameters,   the  partial  derivatives  with 

jF2 
respect  to   the   parameters  t-tt— |      =  I,   the  identity 

DXa  xa=x0 
matrix  and  thus  the  above  equation  reduces  to  the  form 


(A/PiAi  +  P2)X  +  (A1TP1L1  +  P2L2)  =  0 


This  has  the  solution 


X  =  -(A^PxAx  +  Pa)'1  +  (AjPiLj  +  P2L2) 


The  variance  covariance  matrix  of  the  parameter  estimates  is 
obtained  by  normal  error  propagation  methods  and  is  given  by 

\    =  (VPA  +  Pa)'1 

with   the  a  posteriori  variance   of   unit  weight  by 

m  rp 

A         V    P   V    +V    P   V 

*  2  .    1   11   v2  l  2  2 
°o  —  

nl  +  n2  "  u 


in  which  nx  equals  the  number  of  [signal  time  interval] 
observations,   n2  equals  the  number  of  a  priori     parameter 
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observations,  and  u  equals  the  number  of  parameters  in  the 
adjustment. . . . 

If  the  assumptions  going  into  the  adjustment  are  good, 

then  the  a  posteriori    variance  should  be  close  to  1.0.    This 

statistic  is  an  indicator  of  the  quality  of  the  adjustment. 

The  least  squares  adjustment  technique  is  an  excellent 
tool  for  determining  solutions  of  unknown  variables.  However, 
certain  limitations  do  apply.  The  process  assumes  that  all 
systematic  errors  have  been  removed  or  accounted  for,  and  that 
all  remaining  errors  are  randomly  distributed.  Systematic 
errors  that  still  unaccountably  exist  in  the  observations  will 
bias  the  adjustment  such  that  it  may  attempt  to  distribute  the 
errors  to  all  the  observations  and  shift  the  parameters. 

Other  factors  which  may  degrade  the  adjustment  are: 

*  A  poor  physical  model  such  as  one  that  ignores  scale 
error. 

*  The  incorrect  weighting  of  observations. 

*  Small  residuals  which  may  be  a  result  of  poor  network 
geometry  or  insufficient  observations  rather  than  good 
observations . 

For  additional  information  on  adjustment  computations,  see 
Uotila  (1986)  for  derivations  of  appropriate  expressions  for 
least  squares  adjustments  which  use  a  variety  of  different 
systems  or  groups  of  observations  with  their  associated 
constraints. 

Appendix  A  gives  greater  detail  on  the  least  squares 
adjustment  process  for  this  Pegasus  data  analysis. 
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V.   COMPUTATIONAL  PROCEDURE 

The  data  set  used  for  most  of  the  results  below  was  the 
downcast  data  from  Pegasus  Drop  #157  on  August  3,  1989. 
Occasional  comparison  studies  were  also  made  with  the  #157 
upcast  data,  and  that  of  #158  which  occurred  in  the  same 
vicinity  later  the  same  day. 

As  stated  in  Chapter  II,  Pegasus  data  is  recorded  every  16 
seconds  as  the  instrument  descends  and  ascends  during  a 
deployment.  For  this  study,  the  portion  of  the  depth  profile 
of  Drop  #157  from  about  16  m  to  2200  m  was  used,  providing  306 
Pegasus  records.  This  carried  the  instrument  from  surface 
waters  down  through  the  plane  of  the  transponders. 

It  needs  to  be  stressed  at  the  outset  that  the 
computational  procedures  adopted  in  this  thesis  should  not  be 
considered  as  optimum,  but  rather  were  developed  as  processing 
proceeded  in  order  to  overcome  difficulties  encountered  with 
the  specific  data  set  associated  with  the  CIO  array.  If  the 
problems  ultimately  discovered  with  the  data  set  had  been 
known  at  the  beginning,  then  some  procedural  points  in  the 
data  processing  would  have  been  slightly  revised.  This  aside, 
the  purpose  of  the  study  was  to  demonstrate  the  capabilities 
of  least  squares  estimation  procedures  to  resolve  both  unknown 
transponder  positions  and  Pegasus  velocity  components.  As 
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will  be  seen  in  the  following  discussions,  this  was  more  than 
adequately  achieved. 

Because  of  the  poor  initial  positions  of  the  transponders, 
it  was  decided  to  select  a  small  data  set  of  74  Pegasus 
records  from  the  total  of  306  collected  during  drop  #157  and 
use  these  to  help  provide  improved  transponder  coordinates. 
This  data  set  was  determined  by  dividing  the  full  depth  of  the 
Pegasus  drop  into  36  intervals  and  taking  a  pair  of 
consecutive  records  in  each  interval.    Using  the  a  priori 

standard  deviations  derived  in  Chapter  III,  an  initial 
solution  for  the  74  Pegasus  positions  and  the  four 
transponders  was  obtained.  The  very  small  positional  standard 
deviation  associated  with  Tl  essentially  served  to  hold  this 
transponder  fixed  in  the  resulting  solution. 

As  described  in  Chapter  IV,  {A-^'P1Al  +  P2)_1  is  the  variance- 
covariance  matrix  of  the  adjusted  parameters.  The  upper 
tridiagonal  portion  of  this  matrix  is  stored  in  vector  form 
and  can  be  accessed  to  give  information  on  the  variances  and 
covariances  of  the  adjusted  parameters  as  shown  in  Figure  13. 

For  example,  the  portion  of  this  figure  which  describes 

Pegasus  Position  1  gives  the  variances  of  the  x,  y,  and  z 

coordinates,  oPxl2,  o^2,  and  oPzl2,  respectively.    The  other 

three  elements  in  this  upper  tridiagonal,  oP  P   ,  oP  P   ,  and 

xl  yl    xl  zl 

oP   P   ,  provide  the  covariances  between  the  x,  y,  and  z 

yl  zl 

coordinates  of  the  first  Pegasus  position.  Similarly,  the 
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Figure  13.   Variance-Covariance  Matrix 
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covariances  between  the  coordinates  of  different  positions  are 
provided  in  their  appropriate  locations  as  shown.  The 
variances  and  covariances  of  the  transponder  coordinates  are 
located  in  the  lower  right  hand  corner,  and  among  themselves 
comprise  an  upper  tridiagonal  12  x  12  submatrix. 

The  adjustment  solution  is  also  used  to  compute  observa- 
tional residuals  which  in  turn  are  used  to  compute  an  a  posteriori 

variance  of  unit  weight.   For  the  74  record  run  this  a  posteriori 

variance  turned  out  to  be  0.1196,  a  very  small  number  compared 
to  that  which  should  ideally  have  been  close  to  unity.  This 
a  posteriori     variance  when  multiplied   through   the  variance 

covariance  matrix  enables  the  matrix  to  be  scaled  in  order 
that  it  provide  supposedly  unbiased  estimates  of  the  parameter 
variances  and  covariances.  This  was  done  with  this  data  set 
and  the  resulting  transponder  coordinates  with  their  newly 
estimated  standard  deviations  (of  the  order  of  10-20  m)  taken 
and  used  as  a  priori   information  in  the  final  solution  in  which 

all  3  06  records  were  used  simultaneously.  It  was  felt  that 
this  procedure  would  enable  the  transponder  coordinates  with 
their  associated  accuracy  estimates  to  more  closely  represent 
the  type  of  situation  usually  found  in  a  good  transponder 
network. 

In  retrospect,  however,  it  appears  that  it  may  have  been 
more  appropriate  to  have  run  the  74  point  data  set  with  oTime  = 
0.0005  s  rather  than  the  0.002  s  actually  used.   When,  toward 
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the  very  end  of  this  study,  a  oTime  of  0.0005  s  was  allocated  to 
the  measured  time  delays  in  the  74  point  data  set,  the  a  posteriori 

variance  of  unit  weight  became  0.96,  and  the  resulting 
standard  deviations  on  the  transponder  coordinates 
approximately  40  m.  The  coordinate  solutions  for  the 
transponders  did  not  change  by  more  than  5  m  in  position 
although  the  computed  depth  of  T4  did  increase  by  a  further 
10  m. 

Although  this  second  solution  is  to  be  preferred  over  the 
first,  it  was  felt  that  the  work  and  time  involved  in  altering 
all  the  results  and  documentation  already  completed  was  not 
justified,  given  the  minimal  impact  that  it  would  have  on  the 
final  results.  In  fact,  it  was  clear  that  it  would  not  change 
any  of  the  conclusions  resulting  from  this  study. 

Ideally,  when  dealing  with  data  from  a  number  of  Pegasus 
drops  on  a  single  transponder  network,  it  would  be  best  to  use 
a  small,  but  representative  (in  depth)  data  set  from  each 
drop,  and  merge  these  together  in  a  single  adjustment  to 
provide  an  optimum  set  of  transponder  coordinates  for  that 
array  which  could,  where  appropriate,  be  held  fixed  in 
subsequent  simultaneous  processing  of  complete  drops. 
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VI.   RESULTS  OF  STATION  CIO  ADJUSTMENT 

A.   TRANSPONDER  COORDINATES 

With  improved  transponder  coordinates  and  their  associated 
standard  deviations  determined  through  the  74  record  procedure 
described  in  the  last  chapter,  the  entire  306  record  data  set 
was  then  run.  The  results  from  this  full  analysis  are 
discussed  below. 

The  final  transponder  positions  were  moved  (in  total)  from 
their  original  Loran  estimates  as  shown  in  Figure  14  and  Table 
2  below. 
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Figure  14.   Transponder  Movement 
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TABLE  2 
MOVEMENT  OF  TRANSPONDER  POSITIONS 


X 

Y 

Z 

0 

0 

4.5 

-20.3 

-41.9 

-6.2 

32.2 

116.0 

-7.1 

-90.8 

149.3 

-6.5 

Transponder 

Tl 

T2 

T3 

T4 


It  must  be  stressed,  however,  that  this  data  proved  to  be 
only  internally  consistent.  Using  the  same  transponder 
positions  (as  determined  by  the  74  Pegasus  position  data  set) 
to  adjust  data  from  another  drop  #158  (which  used  only  13 
Pegasus  positions) ,  physically  nearby  and  later  the  same  day, 
produced  somewhat  different  transponder  coordinates  although 
it  converged  within  itself.  The  resulting  Table  3  is  shown 
below. 

It  is  suspected  that  the  reason  for  this  inconsistency 
lies  both  in  the  very  low  degrees  of  freedom  in  the  adjustment 
of  drop  #158  (leading  to  a  statistically  weak  solution)  and  in 
the  overall  weak  a  priori   positions  of  the  transponders.   In  a 

well-surveyed  four  transponder  array,  this  problem  would  not 
exist. 

It  appears  certain  that  the  location  of  Pegasus  within 
this  weak  transponder  array  has  an  effect  on  the  final 
adjusted  transponder  positions.  With  the  origin  held  fixed  in 
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TABLE  3 
TRANSPONDER  COORDINATE  SOLUTIONS 

Trans- 


ponder 

X 

Y 

Z 

Tl 

0 

0 

1882.86 

Drop  #157 

T2 
T3 

-1387. 53 
-1769.38 

1530.88 
-66.52 

1987. 15 
2229. 11 

T4 

-540. 12 

2049.71 

1893.48 

Tl 

0 

0 

1879. 13 

Drop  #158 

T2 
T3 

-1244.88 
-1779.68 

1502.60 
-181.94 

1991.22 
2232. 52 

T4 

-533.65 

1950.66 

1896.60 

Tl 

0 

0 

-3.73 

Difference 

Between 

#157-#158 

T2 
T3 

-142 .65 
10.3 

28.28 
115.42 

-4.07 
-3.41 

T4 

-6.47 

99.05 

-3.  12 

both  cases,  the  array  tended  to  skew  in  the  direction  of  the 
drop  (see  Figure  15)  .  A  set  of  coordinates  which  would  be 
appropriate  for  all  drops  might  be  obtained  by  making  a  series 
of  Pegasus  drops  out  along  the  edges  of  the  array  near  the 
centers  of  the  baselines,  in  addition  to  the  drop  in  the 
center. 
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Figure  15.   Adjusted  Array  from  #157  and  #158  Drops 

It  is  well  to  note,  however,  how  the  transponder  positions 
have  been  improved  over  those  of  our  originally  assumed 
coordinates,  especially  those  of  T4  whose  position  was  largely 
unknown  (see  Table  4)  . 

The  horizontal  drms  values  were  (when  using  oT  = 0.0005s): 

*  Tl  =  0  (held  fixed) . 

*  T2  =  11.2  m. 

*  T3  =  13.2  m. 

*  T4  =  12.4  m. 

These  drms  values  suggest  an  improved  horizontal  precision 
of  13.5  meters  or  less. 
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TABLE  4 
STANDARD  DEVIATION  OF  TRANSPONDER  POSITIONS 


Trans 

ponder 

X 

Y 

Z 

Tl 

0 

0 

10 

Initial 

T2 

100 

100 

10 

Standard 

Deviations 

T3 

100 

100 

10 

T4 

100 

100 

100 

Tl 

0.  01 

0.01 

1.88 

Standard 

T2 

8.32 

7.48 

1.75 

Deviations 

After 

T3 

6.79 

11.27 

2.03 

Adjustment 

T4 

11.67 

4  .  12 

2.78 

B.   PEGASUS  POSITIONS 

The  precision  of  the  Pegasus  positions  is  determined  from 
the  variance-covariance  matrix  as  described  in  Chapters  IV  and 
V.   Typical  values  at  various  depths  are  shown  in  Table  5. 

If  geometric  studies  are  desired,  these  variance- 
covariance  values  produce  a  tri-axial  error  ellipse  for  each 
Pegasus  position.  For  ease  of  perception,  the  axes  of  these 
ellipses  can  be  converted  to  an  orthogonal  system  through  the 
use  of  eigenvalues  and  eigenvectors.  (For  a  discussion  of 
this  process,  see  Mikhail,  1976.) 


45 


TABLE  5 
PRECISION  OF  PEGASUS  POSITIONS 


Ox 

oY 

Oz 

4.92 

6.32 

0.51 

4.61 

6.18 

0.49 

4.31 

6.00 

0.50 

3.96 

5.80 

0.53 

3.70 

5.62 

0.58 

3.44 

5.56 

0.67 

3.31 

5.52 

0.86 

3.24 

5.16 

1.29 

3.53 

5.06 

1.98 

3.80 

4.76 

1.53 

Depth  (m) 

16 

248 

503 

750 

1000 

1250 

1502 

1750 

2005 

2183 


C.   CURRENT  VELOCITIES 

Pegasus  velocities  were  determined  by  using  the  adjusted 
values  of  the  Pegasus  positions  which  were  determined  by  the 
methods  described  in  Chapter  V. 

It  is  well  to  note  here  the  similarities  and  differences 
between  the  Pegasus  positions  as  determined  by  this 
adjustment,  and  those  determined  by  the  method  currently  in 
use  by  the  NPS  Oceanography  Department. 

The  initial  Pegasus  positions  used  to  start  this  adjust- 
ment had  been  derived  by  combining  three  observation  equations 
into  two,  and  then  solving  the  two  equations  for  the  unknown 
x  and  y  of  the  Pegasus  position.  (See  procedure  explained  in 
Appendix  A.)    These  positions  were  then  refined  by  the 
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adjustment  process  using  four  transponders  where  Pegasus  depth 
was  constrained  by  pressure. 

The  present  method  used  by  NPS  as  described  in  Chapter  II, 
uses  only  two  observation  eguations  to  solve  for  the  positions 
(i.e.,  using  only  two  transponders),  without  further 
adjustment  and  without  the  pressure/depth  constraint. 
Velocities  derived  from  these  two  sets  of  positions  are 
compared  in  Figures  16  and  17.  (Both  sets  of  velocities  have 
been  filtered  by  a  seven  point  running  average.)  Figure  16 
plots  the  U  components  of  adjustment  derived  velocities 
against  depths,  those  derived  by  using  the  present  NPS  method, 
and  finally  the  differences  between  the  two  approaches. 
Similarly,  Figure  17  plots  comparable  information  for  the  V 
components. 

The  profiles  of  the  two  methods  exhibit  very  similar 
configurations  through  much  of  the  range,  with  notable 
exceptions  occurring  at  the  surface,  at  mid-range  (see 
discussion  of  residuals  in  the  section  on  depth/pressure 
analysis  which  follows) ,  and  at  the  bottom  of  the  cast. 

The  differences  in  surface  velocities  between  the  two 
methods  may  be  a  result  of  using  different  sound  velocity 
profiles.  The  large  divergences  at  depth  where  the  adjusted 
speeds  peak  at  about  22  cm/s  and  those  of  NPS  at  111  cm/s,  are 
almost  certainly  due  to  a  loss  of  precision  in  the  two- 
transponder  solution,   as  well  as  inconsistencies  in  the 
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original  pressure  measurements.  Table  6  provides  a  sample  of 
the  velocities  derived  by  the  two  approaches,  and  shows  the 
differences  between  them.  These  differences  begin  to  increase 
markedly  at  about  2100  m,  indicating  the  substantial  role  the 
adjusted  transponder  network  and  the  inclusion  of  a  pressure 
observation  play  in  the  solution. 


TABLE 

6 

VELOCITY  COMPARISONS 

Adjustment 

Velocity 

Differ- 

Differ- 

derived 

NPS 

Differ- 

ence 

ence 

Depth 

velocity* 

Velocity 

rence 

in  U 

in  V 

(Approx. ) 

(cm/s) 

(cm/s) 

(cm/s) 

(cm/s) 

(cm/s) 

(m) 

m 

(2) 

m-m 

40 

14.23 

13.24 

-0.99 

+  6.30 

+  5.49 

112 

18.89 

25.41 

+  6.52 

+7.09 

+  6.29 

660 

2.10 

1.65 

-0.45 

-0.35 

+  0.29 

1280 

10.01 

10.74 

+  0.74 

-1.01 

-0.77 

1380 

1.08 

9.45 

+  8.37 

-7.06 

-7.74 

1820 

10.34 

12.66 

+  2.33 

+  1.93 

+  1.30 

2090 

17.32 

19.49 

+  2.16 

+  0.81 

+  2.20 

2112 

19.58 

35.84 

+16.27 

+11.61 

+11.44 

2145 

23  .07 

68.05 

+44.98 

+33.89 

+30.10 

2175 

16.36 

110.75 

+94.39 

+71.56 

+62.95 

2190 

16.36 

49.76 

+33.40 

+25.57 

+22.62 

*Adjustment  using  oL  =  0.002,  o   =  2. 


Seven  point  filter  used  for  both  sets  of  velocities. 
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D.   PEGASUS  VELOCITY  ERRORS 

Estimated  velocity  errors  were  obtained  through  the 
propagation  of  positional  errors  found  in  the  variance- 
covariance  matrix  of  the  final  run,  yielding  ou,  ov,  and  ow,  the 
errors  in  the  x,  y,  and  z  directions,  respectively.  This 
propagation  of  errors  proceeds  as  follows. 

The  desired  velocity  components  are  given  by  the 
expression: 


V  = 


U 
V 
W 


(x2-x1)/At 

(y2-yi)/At 

(z2-z1)/At 


where  (x^yuzj  and  (x2/y2/z2)  are  the  positions  of  Pegasus  at 
respective  16  second  intervals  (i.e.,  At  =  16  s) . 
The  variance-covariance  of  V  is  given  by 


2V  = 


uv 
2 


vw 
.  2 


—  geapg 


where 


a/ax,   a/ay-,   a/az,     a/ax2  a/ay2  a/az 


g  = 


l 

At 


-1 
0 
0 


0 

-1 

0 


0 
0 

-1 


1 

0 
0 


0  0 

1  0 
0  1 
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and 


AP 


ZP. 


YP    P 


TP    P 


EP. 


In  this  last  expression  EP,  and  EP^  are  the  variance- 
covariance  matrices  for  Pegasus  positions  1  and  2, 
respectively,  while  2p  P  and  2  P^P,  are  their  cross 
covariances.  These  are  obtained  from  the  full  adjustment 
variance-covariance  matrix  described  in  Chapter  V.  (Appendix 
B  contains  an  algorithm  for  computing  the  variance-covariance 
matrix  of  the  current  velocities  from  data  obtained  from  the 
subroutine  MAT  in  the  Least  Squares  Adjustment  program.) 

These  values  provide  information  on  the  precision  of 
Pegasus  velocity  at  each  position,  and  change  according  to  the 
geometry  between  Pegasus  and  the  transponders.  Table  7 
provides  velocities  and  their  standard  deviations  at  various 
depths  as  computed  by  the  adjustment.  The  horizontal  velocity 
errors  are  greatest  near  the  surface,  and  improve  as  Pegasus 
descends  toward  the  transponders  where,  in  the  horizontal 
sense,  the  geometry  becomes  tighter.  For  the  vertical 
velocity,  the  reverse  is  true. 

These  velocity  errors  were  disappointing.  At  the  surface, 
a  velocity  in  the  X  direction  of  -17.8  cm/s  had  a  standard 
error  of  ±10.7  cm/s.   In  the  Y  direction,  a  velocity  of 
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TABLE  7 
VELOCITIES  AND  STANDARD  DEVIATIONS 


Av. 

Depth 

(m) 

U 

(m/s) 

V 

(m/s) 

W 

(m/s) 

ou 

(m/s) 

°v 

(m/s) 

ow  (m/s) 

20 

-0. 

1784 

0. 

,0169 

0. 

,4347 

0. 

1073 

0. 

,0822 

0.0403 

230 

-0. 

,0886 

0. 

,0789 

0. 

,4566 

0. 

0988 

0. 

,0759 

0.0408 

740 

0. 

,0236 

0. 

,0204 

0. 

,4511 

0. 

0795 

0. 

,0621 

0.045 

1233 

-0. 

,0729 

0. 

,0011 

0. 

,4141 

0. 

0643 

0. 

,0511 

0.0573 

1739 

0. 

,1178 

-0. 

,0238 

0. 

,4734 

0. 

,0578 

0. 

,0439 

0.1107 

2180 

0. 

,1989 

0. 

,0681 

0. 

.3379 

0. 

0569 

0. 

,0454 

0. 1339 

1.7  cm/s  had  a  standard  error  of  ±8.2  cm/s.  At  the  bottom 
results  were  somewhat  improved  with  u  =  19.9  cm/s  and  ou  =  5.7 
cm/s,  v  =  6.8  cm/s  and  ov  =  4.5  cm/s,  but  not  substantially  so. 

These  velocities  and  error  estimates  were  obtained  by 
looking  at  only  individual  pairs  of  positions.  The  NPS 
Oceanography  Department  computes  velocities  by  using  a  seven- 
point  running  average  for  each  Pegasus  recorded  depth.  This 
same  filtering  technique  used  with  the  adjusted  data  would 
improve  precision  by  1/^T  or  0.38,  if  positions  were 
uncorrelated.  However,  these  positions  are  correlated,  so  the 
expected  improvement  is  somewhat  smaller. 

There  are  two  ways  the  precision  of  the  unfiltered 

velocities  could  be  improved: 

*  Use  longer  time  intervals,  so  that  as  At  becomes  larger 
in  the  error  propagation  equation  discussed  above,  the 
variance-covariance  matrix  2  V  diminishes. 
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*  Average  the  16-second  velocities  (in  the  same  manner  as 
NPS)  by  using  a  seven-point  running  average  which  would 
improve  the  precision  if  the  velocities  were 
uncorrelated. 

However,  both  of  these  techniques,  while  gaining  precision, 

sacrifice  vertical  resolution. 

There  are  two  important  aspects  to  these  formal  velocity 
errors.  In  the  first  instance  the  signal  time  intervals  are 
only  measured  to  0.0001  s.  At  a  conventional  sound  velocity 
of  1500  m/s  this  is  equivalent  to  15  cm  in  the  derived  range. 
This  in  turn  propagates  into  a  velocity  error  of  approximately 
1  cm/sec  between  consecutive  Pegasus  positions  if  the 
transponder  positions  are  assumed  to  be  without  error. 

In  the  second  instance,  the  formal  errors  in  the 
transponder  positions  propagate  directly  into  errors  in  the 
Pegasus  positions  and  thence  into  the  velocity  components. 
These  latter  errors  have  by  far  the  most  significant  influence 
on  the  derived  standard  deviations  of  the  velocity  components 
for  Pegasus.  This  in  turn  provides  an  added  emphasis  on  the 
need  for  a  strong  transponder  survey. 

E.   DEPTH/ PRESSURE  ANALYSIS 

To  determine  at  what  depth  the  pressure/depth  measurement 
becomes  critical  in  solving  for  current  velocities,  a  separate 
computer  run  was  made  with  the  standard  deviation  of  the 
Pegasus  Z  coordinate  changed  from  2.00  to  50.00  m.  The 
results  are  graphically  illustrated  in  Figure  18. 
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Figure  18.   Velocities 


Here,  the  U  and  V  velocity  components  derived  from  the 

adjustment  using  oP  =  2  m  are  plotted  on  the  left;  speed 

z 

differences  between  the  U  and  V  components  of  the  two  sets  of 

adjustments  using  oP  =2,  and  oP  =  50,  respectively,  are 

z  z 

illustrated  on  the  right.  The  differences  are  minimal  down 
through  the  water  column  until  about  1700  m  as  Pegasus 
approaches  the  plane  of  the  transponder.  DU  reaches  a  maximum 
positive  difference  of  1.16  cm/s  at  2008  m,  and  subsequent 
maximum   negative   difference   of   -3.53  cm/s   at   2135  m. 
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Respective  values  for  DV  are  0.49  cm/s  at  2001  m  and 
-2.26  cm/s  at  2135  m. 

An  analysis  comparing  the  standard  deviations  of  Pegasus 
depth  positions  shows  a  large  maximum  of  36  m  at  a  depth  of 
2000  m  as  derived  from  the  free-floating  P2  adjustment, 
compared  with  a  maximum  of  2  m  at  a  depth  of  around  194  0  m  for 
tightly-held  P2.   The  latter  reflects  the  fact  that  the  a  priori 

precision  estimate  for  pressure  observations  was  itself  2  m, 
and  that  travel  time  measurement  geometry  gives  no  additional 
information  on  the  adjusted  z  value. 

Table  8  shows  observational  residuals  after  both  of  the 
adjustments  and  compares  the  two  sets  of  residuals  at  similar 
depths.  (These  figures  compare  runs  made  with  oTime  =  0.002  s, 
rather  than  the  final  0.0005  s.) 

At  the  surface  both  sets  of  residuals  are  small.  Where 
the  pressure/depth  relationship  is  tightly  constrained, 
residuals  were  highest  at  the  bottom  of  the  cast  where  Pegasus 
moved  through  the  plane  of  the  transponders,  i.e.,  between 
about  1870  m  to  2240  m.  Here  the  resolution  of  depth  becomes 
less  certain  due  to  the  poorer  geometry  of  the  Pegasus- 
Transponder  array.  These  large  residuals  at  depth  indicate 
there  may  be  a  bias  in  the  pressure  measurements,  thus 
suggesting  the  presence  of  a  systematic  error  in  the  pressure 
sensor  (perhaps  due  to  a  hysteresis  in  temperature  readings) 
that  is  unaccounted  for  in  the  adjustment  model. 
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TABLE  8 
OBSERVATIONAL  RESIDUALS 


Tl 


T2 


T3 


T4 


Surface 

2 

-0. 0001 

-0. 0002 

•  0.0001 

0.0001 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0. 0001 

-0.0001 

-0.0001 

50 

0.0000 

-0.0002 

0. 0001 

0.0001 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0. 0001 

0.0000 

0. 0000 

1400  m 

2 

-0.0004 

-0.0012 

0.0006 

0.0011 

(approx. ) 

-0.0005 

-0.0013 

0.0006 

0.0012 

-0.0005 

-0.0012 

0.0006 

0.0011 

50 

-0.0004 

-0.0011 

0.0006 

0.0010 

-0.0004 

-0.0012 

0.0006 

0.0012 

-0.0004 

-0. 0013 

0.0007 

0.0012 

Bottom 

2 

0.0014 

0.0011 

0.0007 

0.0003 

0.0016 

0.0014 

0.0008 

0.0004 

0.0019 

0.0016 

0.0008 

0.0003 

50 

0.0001 

0.0001 

0.0000 

0.0000 

0.0002 

0.0002 

0.0000 

-0.0001 

0.0002 

0. 0003 

-0.0001 

-0.0002 

Residuals  at  the  bottom  of  the  cast  from  the  adjustment 
where  depth  was  allowed  a  very  loose  constraint  are  uniformly 
low.  The  differences  between  these  two  sets  of  residuals  at 
the  bottom,   as  well  as  differences  in  velocities  and  an 

increase  in  oP   at  this  depth  as  discussed  above,  all  tend  to 

z 

corroborate   the   suggestion   of   a   bias   in   the   pressure 
measurements. 
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Residuals  also  increased  in  both  cases  for  a  few  positions 
around  the  depth  of  1400  meters.  This  appeared  to  occur  where 
the  trajectory  of  the  instrument  reached  its  most  westerly 
position.  Reasons  for  this  are  unclear,  although  a  possible 
explanation  may  be  related  to  the  fact  that  this  depth  and 
position  approximate  the  location  of  the  shelf  edge  of  the 
canyon. 

It  is  also  possible  that  one  or  more  of  the  travel  time 
observations  in  this  part  of  the  drop  had  large  errors.  This 
interpretation  is  further  suggested  by  the  anomalous  current 
shear  computed  at  this  depth  from  the  original  data  (see 
Figures  16  and  17)  . 

It  should  additionally  be  noted  that  these  large  residuals 
at  depth  and  at  mid-range  seem  to  display  a  systematic 
structure  wherein  their  sign  is  consistently  the  same.  If  the 
errors  were  truly  random,  the  positive  and  negative  signs  of 
the  residuals  would  be  distributed  more  or  less  evenly.  This 
situation  is  another  indication  of  some  systematic  error  which 
has  not  been  accounted  for  in  the  adjustment  model. 
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VII.   CONCLUSIONS  AND  RECOMMENDATIONS 

It  is  unfortunate  that  the  data  set  and  the  transponder 
array  used  in  this  study  lacked  the  quality  desirable  for  de- 
tailed analysis.  Many  of  these  difficulties  were,  unfortu- 
nately, only  discovered  as  data  processing  proceeded.  It  is 
encouraging  to  note,  however,  that  the  least  squares  proce- 
dures described  here  enabled  the  identification  of  problems 
with  the  data  which  would  otherwise  have  passed  undetected. 

From  a  theoretical  standpoint,  the  least  squares  process 
must  provide  better  results  than  the  standard  NPS  positioning 
techniques  because  of  the  fact  that  it  uses  all  the  observa- 
tional data.  It  has  the  added  advantage  of  providing  vital 
statistical  information  on  the  quality  of  the  derived  results. 
Unfortunately,  in  the  case  of  the  CIO  network,  the  uncertainty 
regarding  the  original  positions  of  the  transponders  made  it 
difficult  to  perform  a  comprehensive  comparison  between  the 
least  squares  technique  and  the  standard  NPS  methodology. 

It  is  observed,  however,  that  the  least  squares  method 
obtained  a  solution  for  transponder  positions  that  converged 
to  less  than  one  meter,  with  standard  deviations  that  were 
much  smaller  than  those  with  which  the  adjustment  started. 
Precision  of  the  transponder  coordinates  showed  horizontal 
drms  values  of  less  than  15  meters,  with  standard  deviations 
of  transponder  depth  less  than  three  meters.    While  these 
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figures  may  be  optimistic  (see  Chapter  IV) ,  they  are 
indicative  of  the  improvements  which  can  be  achieved  using 
these  methods. 

Regarding  the  precision  of  the  horizontal  Pegasus 
positions,  the  standard  deviation  of  the  X  values  ranged  from 
4.9  down  to  3.2  m,  reaching  a  minimum  around  a  depth  of  1800  m 
just  at  the  upper  limit  of  the  transponder  plane,  and  then 
increasing  slightly  to  the  bottom  of  the  cast.  Pegasus  Y 
values  showed  standard  deviations  of  6.3  to  4.8  m,  with 
passage  through  the  transponder  plane  not  as  noticeable. 

Standard  deviations  of  the  Pegasus  Z  values,  when  they 
were  tied  tightly  to  pressure,  ranged  from  1.4  m  near  the 
surface  to  a  maximum  of  2.0  m  at  a  depth  of  194  3  m.  When  the 
depth  was  only  loosely  constrained,  standard  deviations  varied 
from  2.3  m  near  the  surface,  to  8.6  around  depths  of  1800  m 
just  prior  to  entering  the  plane  of  the  transponders,  and 
reached  a  maximum  of  3  6.5  at  a  depth  of  2  000  m.  From  there 
they  steadily  diminished  back  to  10.4  at  the  final  depth  of 
2165  m. 

Comparison  of  the  results  of  the  NPS  method  and  the  least 
squares  method  shows  a  considerable  difference  in  velocities 
at  depths  where  Pegasus  passes  through  the  plane  of  the 
transponders.  This  difference  at  depth  occurs  also  in 
velocity  comparisons  between  adjustments  made  holding  depth 
constrained  with  a  standard  deviation  of  2.0  m,  and  allowing 
it  to  float  more  freely  with  a  standard  deviation  of  50.0  m. 
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This  leads  to  a  suspicion  that  there  may  be  a  systematic  error 
in  the  pressure  reading  that  has  not  been  accounted  for  in  the 
model,  perhaps  due  to  a  hysteresis  in  temperature  readings. 
Accuracy  of  pressure  observations  becomes  critical  at  depths 
below  1700  m. 

It  was  disappointing  to  note  that  the  standard  deviations 
on  the  computed  Pegasus  velocities  (when  using  At  =  16  sec) 
were  often  of  a  similar  magnitude  as  the  velocities  themselves 
(i.e.,  5-10  cm/s) .  However,  with  resolution  of  signal  travel 
time  possible  only  to  the  nearest  0.0001  s,  resolution  of 
range  becomes  ±15  cm.  This  in  turn  propagates  into  velocity 
errors  in  the  order  of  1-2  cm/sec.  Therefore,  it  is  unrealis- 
tic to  look  for  much  greater  precision  than  that. 

It  is  recommended  that: 

*  NPS  refine  the  existing  least  sguares  technigues  and  the 
adjustment  software  to  facilitate  its  use  on  a  regular 
basis,  for  production  operations. 

*  Four  transponders  be  used  in  each  network.  This  will 
both  strengthen  the  solutions  for  the  Pegasus  velocities 
and  provide  a  reasonable  measure  of  redundancy  for  the 
adjustment  procedure. 

*  More  acceptable  positioning  of  transponders  be  under- 
taken, within  the  time  constraints  involved,  including 
observation  of  all  baselines.  Close  attention  should  be 
paid  to  both  length  and  azimuth. 

*  That  a  small,  but  well-distributed  data  set  of  Pegasus 
records  be  used  from  each  drop  to  assist  in  providing  a 
unique  set  of  transponder  positions  for  each  network  and 
that  these  positions  be  held  fixed  in  the  subsequent 
adjustment  of  the  full  set  of  records  from  each  drop. 

*  Where  practical,  pressure  information  always  be 
collected.   It  both  adds  to  the  overall  redundancy  in  the 
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network  and  enables  additional  parameters  to  be 
introduced  if  necessary. 

A  very  well-controlled  set  of  experimental  data  be 
collected  on  a  well-positioned  five  transponder  network. 
This  will  enable  clarification  of  the  following  issues: 

The  lack  of  consistency  between  the  adjusted  transpon- 
der positions  when  using  different  data  sets  for  the 
same  area. 

The  question  of  pressure/depth  relationships, 
especially  at  depth. 

Possible  causes  for  high  residuals  at  a  given  depth 
(1400  meters  in  the  case  of  station  CIO) . 

-  The  introduction  of  additional  parameters  to  describe 
possible  mis-calibration  of  the  pressure  head  or 
hysteresis  of  temperature  readings. 

-  The  investigation  of  the  actual  motion  of  Pegasus  as 
it  descends  and  ascends. 

The  geometric  impact  of  drops  made  on  the  edges  of  the 
transponder  array  as  opposed  to  those  made  at  the 
center  of  the  array. 
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APPENDIX  A 
LEAST  SQUARES  ADJUSTMENT  PROGRAM 

The  least  squares  adjustment  equation  for  a  combined 
system  of  two  sets  of  observations  is  (Uotila,  1986,  p.  97): 


x    =    -(A/PjAj    +    A2TP2A2)"1(A1TP1L1    +    A2TP2L2) 


where,  in  the  terminology  of  the  program: 

a 

x    =  column  vector  of  corrections  to  the  parameters 

(np,i) ; 

L    =  Observations  (N0BS,1); 

A    =  Jacobian  matrix,  or  the  partial  derivatives  of  each 
observation  with  respect  to  each  parameter 
(NOBS,NP) ; 

A7   =  Transpose  of  A  (NP,NOBS); 

P    =  Weight  matrix  (NOBS, NOBS) ; 

N    =  Maximum  #  of  Pegasus  positions  to  be  determined; 

NT   =  Maximum  #  of  transponders  in  the  network; 

NOBS  =  Maximum  #  of  observations 

NP   =  Maximum  #  of  parameters  allowed  =  (N  x  3)  + 
(NT  x  3)  . 
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A.   FIRST  SET  OF  OBSERVATIONS 
1 .   Lj  Matrix 

The  first  set  of  observations  uses  one-way  travel 
times  from  the  transponders  to  Pegasus  as  the  observations. 
Parameters  are  the  x,  y,  and  z  coordinates  of  each  Pegasus 
position. 

^  =  FA(X) 

Time  =  Distance/Velocity 
Position  1: 


Time!  =  [(XPx-XJ2  +  (YP^Y^2  +  (ZP^ZJ  2]  */V 


Time2  =  [(XPi-X;,)2  +  (YP^Y^2  +  (ZP1-Z2)2]Vv 


Time3  =  [(XPi-Xa)2  +  (YP1-Y3)2  +  (ZP1-Z3)2]Vv 


i 
Time,  =  [(XP^XJ2  +  (YP^YJ2  +  (ZP1-ZJ2]2/V 


Position  N: 

Timei    =  time  intervals  from  transponders  1,4; 

V       =  sound  velocity; 
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XP,YP,ZP   =  Pegasus  position  coordinates; 
X1/Y1,Z1    =  transponder  coordinates; 
#  Observations  =  N  x  4. 

Each  Pegasus  position  therefore  has  four  time-interval 
observations,  one  for  each  transponder.  Thirteen  Pegasus 
positions  will  provide  52  time  observations  and  (13  x  3) 
parameters,  and  therefore  13  degrees  of  freedom  for  this  first 
set  of  observation  equations.  For  each  additional  Pegasus 
position  added  to  the  data  set,  an  additional  degree  of 
freedom  is  gained.   If  a  priori   constraints  are  introduced  for 

the  transponder  coordinates  and  the  transponder  positions 

solved  for  in  the  adjustment  (as  has  been  done  here) ,  then  the 

number  of  observations  and  the  number  of  unknowns  rises  by  12. 

The  adjustment  procedure  requires  an  a  priori   knowledge 

of  the  transponder  positions,  the  depth  of  each  Pegasus 
position  (derived  from  the  pressure  information)  and  the  sound 
velocity  for  each  Pegasus  position.  Estimates  of  the  x  and  y 
coordinates  are  obtained  by  taking  the  three  observation 
equations  for  each  Pegasus  position  that  contain  information 
from  transponders  Tl,  T2 ,  T3  (presumably  the  best  known). 
Subtracting  the  first  equation  from  the  second,  and  the  second 
from  the  third,  will  yield  two  equations  in  two  unknowns,  and 
therefore  a  resulting  first  guess  solution  for  XP  and  YP  at 
that  position. 
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The  Lj  matrix  in  the  adjustment  equation  is  a  matrix 
of  differences  between  observed  time  intervals  and  what  the 
observation  equations  would  yield  usinq  our  first  quess 
parameters. 

■L>1     =     Lx   ca].cuiated     ~     "1   observed 

2 .   A1  Matrix 

The  Ax  matrix  contains  the  derivatives  of  each 
observation  with  reqard  to  each  parameter: 


Ax  Matrix  = 


3x 

a 


If 


Ta  =  [(XP-XJ2  +  (YPj-YJ2  +  (ZPj-Zi)2]  VV  =  D/V 


i  =  1,2  . .  4 
j  =1,2  . .  N 


Then 


5Ti     (XP  -X  ) 
8XPn  = 


1       DV 
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3T 


3ZP 


N 


(ZIW 
DV 


The  derivatives  of  the  transponder  positions  will  be 
the  negatives  of  the  corresponding  Pegasus  position 
derivatives . 


3F. 

i 

3X. 


(xpj-xi) 

DV 


This  will  create  a  matrix  in  the  form: 


Peg. 
Pos  1 


Peg. 
Pos  2 


Peg. 
Pos  N 


3/3p.  d/dp2 


X  ,  X  ,  X 

X  ,  X  ,  X 

X  ,  X  ,  X 

X  ,  X  ,  X 

X  ,  X  ,  X 
X  ,  X  ,  X 
X  ,  X  ,  X 
X  ,  X  ,  X 


3/3T   3/3T„  3/3T-,  3/3T„ 
/    1      2       3      4 


x  ,  x  ,  x 


x  ,  x  ,  x 


x  ,  x  .  x 


x  ,  x  ,  x 


x  ,  x  ,  x 


x  ,  x  ,  x 


x  ,  x  ,  x 


x  ,  x  ,  x 


X  ,  X  ,  X  X  ,  X  ,  X 

X  ,  X  ,  X 

X  ,  X  ,  X 

X  ,  X  ,  X 


X  ,  X  ,  X 


X  ,  X  ,  X 


X  ,  X  ,  X 
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where: 


ti 

=  Time^* 

Pj 

=   XPj,    YP,,    ZPJ? 

Ti 

=   *m  i     *i#    "i i 

i 

=    1,2    . .    4  ; 

J 

=    1,2     . .    N. 

If  submatrices  are  described  by 


aj  = 


X,  X,  X 

X,  X,  X 

X,  X,  X 

X,  X,  X 


x,x,x, 0,0, 0,0, 0,0, 0,0,0 
0,0,0,x,X,X,0, 0,0, 0,0,0 
0,0,0,0,0,0,X,X,X, 0,0,0 
0,0,0,0,0,0,0,0,0,X,X,x 


Then  the  Ax   matrix  becomes 


a:    0     0    0  Sx 

0     a2    0  0  S2 

0     0     a3  0  S3 

aN  SN 


For  each  Pegasus  position  there  will  be  a  4  x  3  block 
for  the  Pegasus  derivatives  and  a  4  x  12  block  for  the 
transponder  derivatives.  All  other  positions  in  the  matrix 
will  be  zero.    Therefore,  to  drastically  reduce  computer 
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memory  and  time  requirements,  the  program  was  developed  to  run 
in  a  sparse  matrix  form  whereby  only  the  non-zero  elements  are 
used.   To  do  this,  the  whole  A  matrix  is  divided  into  parts: 

*  A — which  contains  the  values  for  the  Pegasus  components. 

*  S — which  contains  those  for  the  transponders. 
3  .   P1  Matrix 

Plf    the  weight  matrix,  is  set  up  as  a  diagonal  matrix 

and  stored  in  column  vector  form,  assuming  no  covariances 

between  observations.   If  it  is  desirable  to  add  covariances, 

then  an  error  propagation  subroutine  can  be  inserted  to  fill 

out  the  weight  matrix  and  the  program  altered  accordingly.   In 

our  case  the  chosen  standard  deviation  for  the  time  interval 

was  0.002  s,  later  changed  to  0.0005  s.   It  was  entered  into 

the  program  in  the  error  propagation  section  where  the  code 

reads: 

Do  3  5  I  =  1,NT 
P(l,l)  =  1.0D0/0.0005DO**2 
3  5  Continue 

B.   SECOND  SET  OF  OBSERVATIONS 
1.   L,  Matrix 

The  second  set  of  observations  corresponds  to  various 
parameters  which  have  been  observed,  in  this  case  the  position 
coordinates  of  the  four  transponders,  and  the  position  of  the 
Z  coordinate  of  each  Pegasus  position  which  allows  the  depth/ 
pressure  relationship  to  be  constrained. 
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L2    =    F(x)2 
Xj    =    X, 

Yi    =    Yi 


X,Y,Z    of   each   transponder 


Z«    =    ZA 

ZPi  =  ZPi         Depth  of  Pegasus  as  computed 

i  =  1,N  from  pressure. 

As  described  above,  we  already  have  first  guess 
estimates  of  all  these  observations.  The  L2  matrix  will 
initially  be  zero  since  for  the  first  solution  ,  the  observa- 
tions equal  the  parameters,  i.e.,  Xl  =  Xx.  However,  with 
further  iterations,  this  will  not  be  the  case. 
2  .   A?  Matrix 

The   A2  matrix   contains   the   derivatives   of   each 
observation  with  regard  to  each  parameter: 


A2  =  3F 


3X 
a 


For  example: 


3x- 

7x7 


3X1 

=  0 


3all  else 
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O's: 


This  produces  the  following  diagonal  matrix  of  l's  and 


3  .   P2  Matrix 

P2  is  a  similar  diagonal  matrix  stored  in  column  vector 
form.  It  weights  all  Pegasus  z  positions  to  allow  for 
tightening  or  loosening  the  depth/pressure  relationship.  This 
chosen  standard  deviation,  in  our  case  a  value  of  2.00,  enters 
the  program  through  an  interactive  request  at  the  beginning  of 
the  program  run.   The  standard  deviation  of  each  transponder 
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position  coordinate  is  read  into  the  program  from  the  original 
data  set.   The  V?    matrix  will  appear  as: 


1/    2 


zl 


1/    2 
z2 


1/   2 

/oxl 


1/   2 
/  a  , 


1 


'/„   2 


z4 


C.   MATRIX  OPERATIONS 

Now  all  the  matrices  are  in  place  and  matrix  operations 
can  begin. 

x  =  -(AjP^  +  A2TP2A2)-,(A1TPlL1  +  A2TP2L2) 


72 


Since  the  A2  matrix  is  a  diagonal  of  0*s  and  l's,  the 
second  term  in  the  first  parenthesis  merely  adds  weight  to  the 
appropriate  diagonal  position  in  the  first  term. 

D.  REQUIREMENTS  FOR  RUNNING  THE  PROGRAM 
1.   Initial  Steps 

Before  the  program  is  run,  the  following  steps  must  be 
taken: 

*  Variables  must  be  dimensioned  adequately,  following  the 
definitions  clearly  stated  in  the  preface  to  the  program. 

*  The  value  of  oTime  must  be  defined  (in  our  case  0.0005), 
and/or  the  appropriate  error  propagation  subroutine 
added. 

*  A  data  set  file  must  be  available  for  access  and  set  up 
in  the  following  fashion: 

Transponder  coordinates,  each  with  its  own  standard 
deviation. 

For  each  Pegasus  position:  pressure,  depth,  four 
transponder  one-way  time  intervals,  and  sound 
velocity.   (Sample  data  set  is  shown  below.) 
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Transponder 

Tl 

0.0 

0.0            1878. 

4 

' 

Coordinates 

T2 

-132^.6 

1485.8         1993. 

4 

T3 

-laoi .o 

-1U2.5         2236. 

i 

T4 

-449.3 

00.0 

1900.4         1900.0 
1         00.01            10.00 

Transponder 

Tl 

Standard 

T2 

50.00         50.00            10.00 

Deviations 

T3 

50.00         50.00            10.00 

T4 

50.00         50.00          100.00 
30.3           30.6       1.1K5       1.5486 

1.6277 

1  .6803 

1485.310 

Pegasus 

1 

Positions 

2 

2  i  3  .  4 

211.7       1.3255 

1  .4245 

1  .5191 

1  .5704 

14a5.  100 

3 

393  .0 

394.7       1.2356 

1  .308  3 

1 .4105 

1  .4712 

1485. 140 

A 

581  .5 

576.4       1.1414 

1  .2047 

1 .3062 

1  .3844 

1485.250 

B 

7  6', .  7 

757.7       1.0457 

1.1140 

1 .2043 

1.3106 

1485.650 

. 

94  4  .0 

934.9      0.9592 

1  .0282 

1  .  1  121 

1  .2403 

1486.200 

1125.8 

1114.5      0.8941 

0.93UI 

1.01/9 

1 . 1 754 

1486.820 

• 

1206. 1 

1292.4      0.8405 

0.8o31 

0 . 9243 

1 . 1302 

14b 7.  570 

K89.7 

1472.5      0.7990 

0.7  944 

0.84  72 

I .0851 

1483  .440 

1667.0 

1648.1       0.7620 

0.7455 

0.7922 

1 .0491 

1489.370 

13^6.2 

1824.5      0.7341 

0.7229 

0. 7  5  34 

1 .0306 

1490.400 

202  1.9 

1997.4      0.7337 

0.6973 

0.7584 

0.9940 

1491  .510 

2200.  7 

2173.1       0.7415 

0.6950 

0.7910 

0.9671 

1492.760 

Press- 

Depth 

Ti 

T2 

Ji 

TA 

Sound 

ure 

Velocity 

Time  Intervals 

» 

*  FILEDEFS  for  input  and  output  files  must  be  in  place, 
either  typed  in  before  the  program  starts  or  available 
through  a  program  exec.  The  program  will  run  in  either 
WF77  or  FORTVS2  fortran. 

-  FILEDEF  01  DISK  fn  ft  fm  (i.e.,  Test  Data  A). 

-  FILEDEF  02  DISK(recfm  vb  lrecl  132  blksize  134) . 

(For  large  data  sets,  use  FILEDEF  02  DISK(recfm  fb 
lrecl  132  blksize  13200) . 

*  Be  sure  enough  memory  is  available.  Two  M  sufficed  for 
running  at  least  up  through  74  Pegasus  positions,  but 
4096  K  were  reguired  for  the  306  position  run. 

The  program  starts  off  by  asking  for  operator  input 

(sample  answers  in  parentheses) ,  reguesting: 
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*  Project  or  run  description. 

*  Number  of  transponders  (4). 

*  Global  standard  deviation  to  be  given  to  the  Pegasus  z 
coordinates  (2.00). 

*  Convergence  limit  for  adjustment  and  number  of  iterations 
(entered  separately)  (1.000),  (4). 

From  then  on  the  program  works  on  its  own,  finally 

producing  whatever  output  has  been  requested  in  the  output 

file. 

2 .   Brief  Program  Outline 

*  Matrices  are  zeroed  out  at  appropriate  times. 

*  Pegasus  data  set  is  read  in. 

*  Subroutine  APPRO  is  called  to  calculate  initial  approxi- 
mate coordinates  of  Pegasus. 

*  The  diagonal  weight  matrix  P1  is  set  up  and  stored  in 
matrix  BIGP. 

*  Program  constants  are  calculated. 

*  Weight  matrix  P2  is  generated. 

*  The  Lj  matrix  is  formed. 

*  The  matrices  A  and  S  are  formed  by  calling  subroutine 
SETUP. 

*  The  normal  matrix  equations  are  formed  block  by  block  by 
calling  subroutine  NORM.  Outputs  include  the  normal 
matrix  (ATPA) ,  in  column  vector  form  called  anorm,  and 
xhat  the  ATPL  vector  called  xhat  (which  subsequently 
becomes  the  solution  vector) . 

*  Subroutine  P2L2  is  called  to  compute  the  P2L2  matrix  and 
add  it  to  the  column  vector.  It  also  increments  the 
diagonal  elements  of  the  normal  matrix  to  include  the 
influence  of  P2 . 

*  The  normal  equations  are  solved  by  calling  a  user 
supplied  routine  that  will  accomplish  the  required  matrix 
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operations  for  producing  the  adjusted  values  with  which 
to  correct  the  parameters.  The  IMSL  routine  LINV3P  was 
used  in  this  study.  Utilizing  a  symmetric  storage  mode, 
this  algorithm  replaces  the  ATPA  matrix  by  its  inverse 
which  is  the  variance-covariance  matrix  of  the  Pegasus 
and  transponder  positions.  This  (ATPA)"1  matrix  can  be 
written  to  file  for  later  access  in  order  to  produce  the 
error  matrix  required  for  current  velocity  analysis  (see 
Chapter  V)  as  used  in  the  program  provided  in  Appendix  B. 

*  The  parameter  corrections  are  tested  to  see  if  iteration 
is  required.  If  so,  the  program  then  updates  the 
parameters  and  goes  through  the  process  again  for  as  many 
iterations  as  are  called  for,  or  until  the  convergence 
limit  has  been  reached. 

*  Residuals  and  the  a  posteriori  variance  of  unit  weight  are 
computed  by  calling  subroutine  RESID.  (Residuals  here 
are,  of  course,  given  in  units  of  seconds  since  they 
relate  to  time  interval  observations.) 

3 .   Program  Outputs 

Outputs  include: 

*  Original  data  set  and  interactive  information  given  at 
the  beginning  of  the  program. 

*  Adjusted  parameters  following  each  interaction. 

*  Residuals  and  a  posteriori   variance  of  unit  weight. 

*  Any  other  information  requested  to  be  written  to  a  file. 

For  this  study  a  piece  of  code,  Subroutine  MAT,  was 
attached  at  the  end  of  the  program  to  pick  out  from  the  anorm 
vector  only  the  data  required  for  current  velocity  analysis. 
This  subroutine  is  called  after  the  write  statement  for  the  a 

posteriori    variance  of  unit  weight  which  uses  format  statement 

#175.  The  anorm  vector  is  the  upper  tridiagonal  of  the 
(ATPA) _1  variance-covariance  matrix  stored  in  column  vector 
form.   (See  the  diagram  in  Chapter  V.) 
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The  position  of  the  variance  of  each  parameter  can  be 
obtained  by: 

N  x  (N  +1) 


where  N  is  the  place  number  of  the  parameter.  For  example, 
suppose  the  variance  for  the  Pegasus  y  coordinate  in  the 
second  record  was  requested.  YP2  is  the  5th  parameter — XPX, 
YPi,  ZP:,  XP2,  YP2/  .  ..,  (5*6)/2  =  15.  The  variance  for  YP2/ 
therefore,  will  occupy  the  15th  place  in  the  anorm  vector. 

XPX    YPj    ZPj    XP2    YP2    ZP2 


1 

2 

4 

7 

11 

16 

Peg. 

Pos. 

1 

3 

5 

8 

12 

17 

6 

9 

10 

13 

14 

18 
19 

Peg. 

Pos.  2 

© 

20 

21 

The  subroutine  Mat  picks  up  only  those  values  in  the 
anorm  vector  that  fill  the  3x3  variance-covariance  matrix  of 
each  Pegasus  parameter,  and  the  3x3  covariance  portion  that 
relates  each  parameter  to  the  next  one  in  line.  A  sample  of 
the  output  follows: 
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Var-covar 
Pos  1 


24.248731431(153141 
18.(32808041824(140 
0. 134 (I  (24 14 3  1 1102 7 (E-01 


18.(3280804I62((I40 

31  .8125418557777184 

-0.S115O57775050IO4O5 


0.1544|(24145111027(E-Ol 
-0.511505777505010405 
0.2(071(185552885254 


Covar 

Pos  i-Pos  2 


22.  7  00  7»2  3*2  788  35  4  4        18. 5171145727731243 

18.54  128  213(0100441         31.03(4  505  1(55(78  10. 

0.8316874(41 15 04B052E- 01  -0. 45 15210(1  184188830 


0.43I2402  3I38440S242E-OI 
■0.444420407711407421 
0.520710122205(5712  IE- 01 


Var-covar 
Pos  2 


24.01117S0884437B42 
18.4018782718410138 
0.72311 74421 4 784888SE-0 I 


18  .4018782718410138 

31.1I0223740S132344 

■0.51501(((SI08(57157 


0.72311  74421  4 78 4888 SE- 01 
-0.51501(4(5108(37157 
0.2(0001015573101355 


Covar 

Pos  2-Pos  3 


22.5(71757113281141 
18.55015075515(8424 
0.  434724100041810147 E- 01 


18.5701158104511001 

31.04141780*2177121 

-0.443411422131407538 


0.485  7  714I8  740151404E-01 
■0.445154214583781211 
0.51 0857 48 7137755141 E- 01 


Var-covar 
Pos  3 


2  3.1448  7  4488  488134  7        18.514  34110428  44541 

18.514  341104284(341        31.10182  341584144(1 

0.5(517123((7(8  3  7824E-01  -0.513118(18135371114 


0.5(51 7  12  5(4 7  48 5  78  24E-01 
-0.515118(18155571114 
0.251115(1077(075011 


The  first  set  of  three  lines  with  three  values  each  is  the 
variance-covariance  matrix  for  Pegasus  position  1,  which  is 
symmetric.  The  second  three  lines  of  three  values  is  the 
covariance  matrix  for  position  1  and  position  2 
(nonsymmetric) .  The  third  set  of  three  lines  is  the  variance- 
covariance  matrix  for  Pegasus  position  2  (symmetric) ,  and  so 
on. 

The  final  30  lines  in  the  output  provide  the  variances 
of  the  transponder  coordinates  and  all  the  covariances  between 
them. 
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c 
c 
c 
c 


C  PROGRAM   PEGASUS 

C 

Q  ***************************************************************** 

C     *  * 

C    *        COMMENTS  GIVING  A  DESCRIPTION  OF  THE  PROGRAM 

C    *  * 

Q  ***************************************************************** 

c 

IMPLICIT  REAL*8  (A-H.O-Z) 
C 

c   @(a@@(a<a(a(a<a(a^ 

C   <a  "  @ 

C  @  SET  MATRIX  VARIABLES  TO  THEIR   MAXIMUM   SIZES.    THE   VARIABLES  @ 

C  @  USED  ARE   AS   FOLLOWS:                                                                                              @ 

C  (§  N             MAX.    //   OF   PEGASUS   POSITIONS   TO   BE   DETERMINED   ( 20    )   @ 

C  @  NT     =   MAX.    //   OF  TRANSPONDERS    IN  THE   NETWORK                      (4)      @ 

C  @  NP     =  MAX.    #   OF   PARAMETERS   ALLOWED  =   (N*3)    +   (NT*3)               @ 

C  @  NOBS=   MAX.    #   OF   OBSERVATIONS   =   N*NT                                                    @ 

C  @  NV      =  MAX.    //   OF  NON-ZERO   VALUES    IN  THE   A   MATRIX  =   6*NT*N  @ 

C  @  NN     =  MAX.    #   OF   ELEMENTS    IN  THE   UPPER  TRIANGULAR   PORTION  @ 

C  (9  OF  THE   NORMAL  MATRIX   =   NP*(NP   +   l)/2                                    @ 

C  @  (? 

C  @  THE   FOLLOWING   MATRICES   MUST   BE   DIMENSIONED  TO  THEIR  MAX.         @ 

C  @  SIZE   PRIOR  TO  THE   PROGRAM   BEING   USED.                                                     @ 

C  @  @ 

C  @  XO(NP),    P(NT.NT),    BIGP(NT,NOBS),LB(NOBS,l),    TRANS(NT),          @ 

C  @  X(NT),Y(NT),Z(NT),A(NT,3),S(NT,NT',0),VALUE(NV))ATP(NT,3),   (3 

C  @  ATPA(3,3),    ATPS(3,NT*3),    ATPL(3,1),    STP(NT*3  ,NT)  ,SV(N)          @ 

C  @  STPS(NT*3,NT*3),    STPL(NT*3 , 1)  ,    Ll(NOBS,l),    L(NT,1),                 @ 

C  @  ICOL(NV),IROW(NV),    SPS(NT*3  ,NT*3)  ,SPL(NT*3, 1)  ,ANORM(NN) ,      @ 

C  @  XHAT(NP),XA(NP),VTP(1,NT),V(NT,1),V1(N0BS),SIG(1,1),          @ 

C  @  SDX(NT))SDY(NT),SDZ(NT),P2(NP))L2(NP),LB2(NP),AUX(2)              @ 

C  Q  @ 

C  @  NOTE:    THE   SAME  TRANSPONDERS   MUST  APPEAR   IN  EVERY  PEGASUS     @ 

C  @  POSITION.    A  RECORD    IN  WHICH  ONE   (OR  MORE)   DROP  OUT     @ 

C  @  IS  NOT  ALLOWED.                                                                                      @ 

C  @  @ 

c 

DIMENSION  XO(  72) ,P(4 ,4) ,BIGP( 4  ,  80) ,TRANS(4) ,X(4) , Y(4) ,Z(4) , 

1  A(4,3),S(4,12),VALUE(480  )  ,  ATP(4 , 3) , ATPA( 3  ,  3)  ,  ATPS( 3 , 12) ,SV(  20)  , 

2  ATPL(3,1),  STP(12,4),  STPS(12,12),  STPL(12,1),  ANORM(  2628), 

3  XHAT(  72),XA(  72) , VTP( 1 ,4) , V(4 , 1) , Vl(  80) ,SIG( 1 , 1) ,  SDX(4), 

4  SDY(4),  SDZ(4),  P2(  72),  AUX(2) 

DOUBLE  PRECISION  LB(  80,1),  Ll(  80,1),  L(4,l),  L2(  72),  LB2(  72) 

CHARACTER*12  PROJ 

INTEGER*2  ICOL(480  ),  IROW(480  ) 

COMMON  SPS(12,12),  SPL(12,1) 


C 
C 

c   <a  @ 

C   @     READ  IN  THE  FIXED  DATA  AND  ANY  ASSOCIATED  VARIANCES    @ 
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c 

@ 

c 

@ 

c 

@ 

c 

@ 

c 

@ 

c 

c 

UNITS  5  AND  6  READ  AND  WRITE  FROM  AND  TO  THE  SCREEN  @ 
FOR  INTERACTIVE  PROCESSING.  UNITS  1  AND  2  READ  AND  @ 
WRITE  FROM  AND  TO  DISK  FILES.  @ 

@ 


WRITE(6,*)  'PROJECT  OR  RUN  DESCRIPTION  (TO  60  CHARACTERS)' 

READ(l.lO)  PROJ 
10  FORMAT(A12) 

WRITE(6,*)  'NUMBER  OF  TRANSPONDERS  (THEIR  COORDINATES  AND  STD.  ', 
1  'DEVIATIONS  TO  BE  READ  FROM  THE  DISK  FILE)' 

READ(1,*)   NT 

READ(1,*)    (X(I),Y(I),Z(I),    1=1, NT) 

READ(1,*)    (SDX(I).SDY(I) ,SDZ(I),    I    =    1,NT) 

WRITE(6,*)    'THE   GLOBAL   STD.    DEVIATION  TO   BE   GIVEN  TO  THE    '  , 
1    ' Z   PEGASUS   COORDINATES' 

READ(1,*)   SDZP 

WRITE(6,*)    'CONVERGENCE   LIMIT  FOR  ADJUSTMENT  AND  //   OF   ITERATIONS' 

READ(1,*)  TOL,  NITER 

WRITE(2,15)  PROJ,NT,(X(I),SDX(I),Y(I),SDY(I),Z(I),SDZ(I),I  =1,NT) 

15  FORMAT('l',///,5X,70("*'),//,10X,A60,//,5X,70('*,),//)10X, 

1  'NUMBER  OF  TRANSPONDERS  =', 12 ,//,' TRANSPONDER  COORDINATES  AND  ', 

2  'THEIR  STANDARD  DEVIATIONS ',// ,4(  10X,  3(F10.  2  ,  2X.F6.  2) ,/),/// ) 
WRITE(2,16)  SDZP 

16  FORMAT( '  ',5X,'THE  GLOBAL  STANDARD  DEVIATION  FOR  THE  PEGASUS  ', 
1  'Z  COORDINATES  =' ,F7. 2) 

WRITE(2,17)  TOL,  NITER 

17  FORMATC  ',//,5X,'THE  CONVERGENCE  LIMIT  FOR  THE  ADJUSTMENT 

1  F7.  3,//,6X,'THE  NUMBER  OF  ITERATIONS  ALLOWED  =' ,  13,///, 

2  50X, ' INPUT  DATA' ,//,5X, 'PRESS. ' ,4X, 'ZP     TIME  DELAY  l',2X, 

3  'TIME  DELAY  2',2X,'TIME  DELAY  3',2X,'TIME  DELAY  4',  5X,' SOUND  ', 

4  'VELOCITY'  ,2X,'APPROX.  PEGASUS  POSITIONS(X, Y, Z) ' ,/) 
C 

C   @@(^(a(a(a(a(a(a(a(<^@^ 

c   @  "  @ 

C   @     READ  IN  THE  PEGASUS  DATA,  POSITION  BY  POSITION         @ 
C   @  @ 

c   @@(a^(a(a(a(a(a(a^ 
c 

C   FIRST  ZERO  OUT  THE  P  MATRIX 
C 

,   DO  20  I  =  l.NT 
DO  20  J  =  1, NT 
P(I,J)  =  0.  0D0 
20  CONTINUE 

N  =  0 
25  READ(1,*,  END  =  45)  PRESS,  ZP , (TRANS( I  )  ,  1  =  1 ,NT) , SV(N+1) 
NNT  =  N*NT 
DO  30  I  =  1, NT 
LB(NNT  +  1,1)  =  TRANS(I) 
30  CONTINUE 
C 

C   CALCULATE  THE  APPROXIMATE  POSITION  OF  PEGASUS 
C 
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CALL  APPR0(TRANS,X,Y,Z,SV(N+1) ,NT, AUX.XP , YP , ZP) 
C 

XO(N*3  +  1)  =  XP 

XO(N*3  +  2)  =  YP 

XO(N*3  +  3)  =  ZP 

WRITE(2  32)  PRESS,ZP,(TRANS(I),I=1,NT),SV(N+1),XP,YP,ZP 
32  FORMAT( '  '  , 3X.F7.  2 , 2X.F7.  2 ,5X ,F7. 4 , 3( 6X.F7. 4) , 13X.F8.  3 , 7X, 
1  3(2X,F7. 1)) 
C 

C    DO  ERROR  PROPAGATION  AND  FORM  THE  P  MATRIX  FOR  THE  FIRST  PEGASUS 
C    POSITION.  STORE  THIS  IN  MATRIX  BIGP.  THE  ERROR  PROPAGATION 
C    SUBROUTINE  MUST  BE  INSERTED  HERE  AND  THE  NEXT  FEW  LINES  OF  CODE 
C   MODIFIED  ACCORDINGLY. 
C 
C  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 

c 

C   FOR  INITIAL  PROGRAM  TESTING,  INSERT  DIAGONAL  VALUES  ONLY 
C 

DO  35  I  =  1, NT 

P(I,I)  =  l.ODO/O.  0005D0**2 
35  CONTINUE 

C 

DO  40  I  =  l.NT 
DO  40  J  =  1, NT 
BIGP(I,NNT  +  J)  =  P(I,J) 
40  CONTINUE 
N  =  N  +  1 
GO  TO  25 

45  CONTINUE 
C 

C    CALCULATE  PROGRAM  CONSTANTS 
C 

NMAX  =  N 

NP  =  N*3  +  (NT*3) 

NOBS  =  N*NT 

NV  =  (6*NT)  *  N 

NN  =  NP*(NP+l)/2 

NT3=  NT* 3 

ITER  =  0 

WRITE(2,46)  NMAX, NOBS, NP 

46  FORMATC'l' ,///, 5 X,' ACTUAL  #  OF  PEGASUS  POSITIONS  =',15,//, 

1  5X,' TOTAL  //  OF  OBS.  ON  PEGASUS   =',15,//, 

2  5X, 'TOTAL  #  OF  PARAMETERS        =',I5) 
C 

c   @<a<a(a(a@(c^@^^ 

c   @  "  @ 

C  @  SET  UP  THE  MATRICES  NEEDED  FOR  THE  ADJUSTMENT,  BEGINNING  @ 

C  @  WITH  THOSE  ASSOCIATED  WITH  THE  "OBSERVED"  PARAMETERS.  @ 

C  @  MOST  OF  THESE  NEED  TO  BE  ZEROED  OUT  AT  THIS  POINT.  @ 

C  @  @ 

C  @  THEN  SET  UP  THE  NORMAL  MATRIX  BLOCK  BY  BLOCK,  EACH  BLOCK  @ 

C  @  CORRESPONDING' TO  A  NEW  PEGASUS  POSITION  @ 

c   <a  @ 

c 
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ZERO  OUT  MATRICES 

DO  50  I  =  l.NP 
LB2(I)  =  0.  ODO 
P2(I)  =  0.  ODO 
L2(I)  =  0.  ODO 
XHAT(I)  =  0.  ODO 
50  CONTINUE 

INSERT  THE  NON-ZERO  ELEMENTS 

J  =  N*3 

DO  55  I  =  3, J, 3 
P2(I)  =  1.0D0/SDZP**2 
LB2(I)=  XO(I) 
55  CONTINUE 


58 


800 


801 


802 


803 


J  =  N*3  +  1 

K  =  1 

DO  58  I 

XO(I) 

P2(I) 

LB2(I)   = 

X0(I+1)  = 

P2(I+1)  = 

LB2(I+1)= 

X0(I+2)  = 

P2(I+2)  = 

LB2(I+2)= 

K  =  K  +  1 

CONTINUE 

WRITE(2,800)(LB(I,1),I=1,NOBS) 

FORMAT( '1'  ,///,50X,*LB  MATRIX' ,// ,8( IX,  10(F10. 7 , 2X)/) , 1X.2F10.  7) 

WRITE(2,801)(XO(I),I=1,NP) 

FORMAT( *0'  ,///,50X,'XO  VECTOR'  , // , 8( IX, 10(F10. 4,2X)/) , 1X.F10.  4) 

WRITE(2,802)(P2(I),I=1,NP) 

FORMAT( '0'  ,///,50X,'P2  VECTOR'  ,// ,8(  IX, 10(F10. 4 ,2X)/) , 1X.F10.  4) 

WRITE(2,803)(LB2(1),I=1,NP) 

FORMAT( '0'  ,///,50X, ' LB2  VECTOR'  , //  ,  8(  IX,  10(F10. 4, 2X)/) , 1X.F10.  4) 


J,NP,3 

X(K) 

1.  0D0/SDX(K)**2 

XO(I) 

Y(K) 

1. 0D0/SDY(K)**2 

X0(I+1) 

Z(K) 

1.  0D0/SDZ(K)**2 

X0(I+2) 


@  @ 

@  THE  ITERATIVE  PROCESS  BEGINS  FROM  HERE.  BEGIN  BY  ZEROING   @ 

@  OUT  THE  NORMAL  MATRIX  AND  THE  TWO  COMMON  BLOCK  MATRICES    @ 

@  @ 

60  DO  62  I  =  1.NT3 

SPL(I.l)  =  0.  ODO 

DO  62  J  =  1.NT3 

SPS(I.J)  =  O.ODO 
62  CONTINUE 

DO  65  I  =  1,NN  ' 

ANORM(I)  =  0.  ODO 
65  CONTINUE 


82 


K  =  1 

J  =  N*3  +  1 
DO  66  I  =  J,NP,3 
X(K)  =  XO(I) 
Y(K)  =  XO(I+l) 
Z(K)  =  X0(I+2) 
K  =  K+l 
66  CONTINUE 
C 

c 

C    BEGIN  THE  FORMATION  OF  THE  A  AND  LI  MATRICES 
C 

DO  100  K  =  l.NMAX 

JJ  =  (K-1)*NT 

DO  68  I  =  1, NT 

TRANS(I)  =  LB(JJ  +  1,1) 

DO  68  J  =  1, NT 

P(I,J)  =  BIGP(I,JJ  +  J) 
68  CONTINUE 

KK  =  (K-l)*3 

XP  =  XO(KX  +  1) 

YP  =  XO(KK  +  2) 

ZP  =  XO(KK  +  3) 
C 
C 

DO  70  I  =  1,NT 

L1(JJ  +  1,1)  =  (DSQRT((XP-X(I))**2  +  (YP-Y(I))**2  +  ( ZP-Z( I ) )**2) 
1  /SV(K))  -  TRANS(I) 

L(I,1)  =  L1(JJ  +  1,1) 
70  CONTINUE 
C 

C    ZERO  OUT  THE  A  AND  S  SUBMATRICES 
C 

DO  80  I  =  1, NT 

DO  75  J  =  1,3 

A(I,J)  =  0.0D0 
75  CONTINUE 

DO  80  J  =  1.NT3 

S(I,J)  =  0. 0D0 
80  CONTINUE 
C 

CALL  SETUP(X,Y)Z,XP,YP,ZP)SV(K),K,NTJNV,NT3,NMAX,A,S,ICOL,IROW, 
1  VALUE) 
C 

CALL  NORM(A,S,P,L,K,NT,NT3,NN,NP,NMAX,ATP,ATPA,ATPS,ATPL,STP, 
1  STPS,STPL,ANORM,XHAT) 
100  CONTINUE 
C 

WRITE(2,804)(L1(I,1),I=1,N0BS) 
804  FORMATC'l'  ,///,50X,'Ll  MATRIX'  , // , 8( IX,  10(F10. 7 ,2X)/) , 1X.2F10.  7) 
C 

C   PLACE  SPS  AND  SPL  INTO  THE  NORMAL  MATRIX  AND  THE  COLUMN  VECTOR 
C   RESPECTIVELY. 
C 

DO  110  I  =  1.NT3 

NCOL  =  (NP  +  1)  -  I 
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11  =  NCOL*(NCOL  +  l)/2  +  I 
DO  110  J  =  I.NT3 

12  =  II  -  J 

ANORM(I2)  =  SPS(NT3-J+1,NT3-I+1) 
110  CONTINUE 

13  =  NMAX*3  +  1 
DO  120  I  =  13, NP 

XHAT( I)  =  SPL( 1-13+1,1) 
120  CONTINUE 
C 

C    SET  UP  THE  L2  MATRIX  AND  ADJUST  THE  DIAGONAL  ELEMENTS  OF  THE 
C    NORMAL  MATRIX  FOR  THE  INFLUENCE  OF  P2 
C 

CALL  P2L2( P2 , LB2 , XO , ANORM , L2 , XHAT , NP , NN ) 
C 

c   @(^(a@(^(a(a<a(a(^^ 

c   <a  @ 

C   <a    SOLVE  THE  NORMAL  EQUATIONS  AND  COMPUTE  THE  RESIDUALS     @ 
C   (a  @ 


C     CALL  A  USER  SUPPLIED  ROUTINE  TO  SOLVE  THE  NORMAL  EQUATIONS. 
C     tSEE  APPENDIX  A,  SECTION  D.) 

DO  135  I  =  l.NP 

XHAT(I)  =  -XHAT(I) 

XA(I)  =  XO(I)  +  XHAT(I) 

XO(I)  =  X*(I) 
135  CONTINUE 

ITER  =  ITER  +  1 

WRITE(2,140) 
140  FORMAT( '0' ,///,5X,' ADJUSTED  PEGASUS  POSITIONS  (WITH  THE  ', 

1  'CORRECTIONS  TO  THE  APPROX.  POSITIONS  IN  PARENTHESES)',//,  18X, 

2  'X' ,18X,'Y' ,18X,'Z') 
J  =  NP  -  NT3 

DO  145  1=1, J, 3 

WRITE(2,142)  XA(I),XHAT(I) ,XA(I+1) ,XHAT( 1+1) ,XA( 1+2) , 
1  XHAT(I+2) 
142  FORMAT( '  '  , 10X.F9. 3 , 2X, ' ( *  ,F6.  1  ,  '  )  '  , 2( 2X,F9.  3 , 2X, ' ( '  ,F6.  1 , ' ) ' ) ) 

145  CONTINUE 
J  =  J  +  1 
WRITE(2,146) 

146  FORMATC  ' ,//,5X,' ADJUSTED  TRANSPONDER  POSITIONS  (WITH*, 
1  'CORRECTIONS  TO  TOE  APPROX.  POSITIONS  IN  PARENTHESES)') 

DO  148  I=J,NP,3 

WRITE( 2 , 142)  XA( I) ,XHAT( I ) ,XA(  1  +  1 ) ,XHAT( 1  +  1 ) ,XA( 1+2)  , 
1  XHAT(I+2) 
148  CONTINUE 
C 

C   TEST  THE  PARAMETER  CORRECTIONS  TO  SEE  IF  ITERATION  IS  REQUIRED 
C 

I FLAG  =  0 
DO  150  I  =1,NP 

IF(DABS(XHAT(I)).GT.  TOL)IFLAG  =    IFLAG   +    1 
150   CONTINUE 
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C    COMMENT  OUT  THE  NEXT  THREE  STATEMENTS  DURING  PROGRAM  DEVELOPMENT 
IF(IFLAG.EQ.  0)  GO  TO  155 

IF(  ITER.  LT.  NITER)  GO  TO  60 
C 

C    COMPUTE  THE  RESIDUALS 
C 

155  CALL  RESIDCVALUE.ICOL.IROW.BIGP.Ll.XHAT.NMAX.NT.NV.NP.NOBS, 
1  P,VTP,V,SIG,V1,SIG0) 
C 

WRITE(2  160) 
160  FORMATC'l' ,///,10X,' OBSERVATIONAL  RESIDUALS  AFTER  THE  ADJUSTMENT', 
1  //,5X, 'TRANSPONDER  1    TRANSPONDER  2     TRANSPONDER  3',4X, 
2 'TRANSPONDER  4'  ) 
DO  170  I  =  1, NOBS, NT 
WRITE(2  165)(V1(I+K-1),  K  =1,NT) 
165  FORMATC   ' , 9X.F8. 4 , 3( 9X.F8. 4) ) 
170  CONTINUE 

WRITE(2  175)  SIGO 
175  FORMAT( '  ',///, 5X ,' THE  A  POSTERIORI  VARIANCE  OF  UNIT  WEIGHT 
1  F8.  4) 
C 

C    *THE  FOLLOWING  WRITE  STATEMENT  WAS  ADDED  BY  M.  HASKELL  TO  LOOK 
C    *AT  THE  VARIANCE-COVARIANCE  MATRIX  FOR  POSITION  COORDINATES. 
C    *THE  CALL  STATEMENT  FOR  SUBROUTINE  MAT  WAS  ADDED  TO  ACCESS  ONLY 
C     "THE  INFORMATION  NEEDED  TO  CALCULATE  THE  VARIANCES  AND  COVARIANCES 
C    *OF  THE  VELOCITIES. 

C  . 

C  WRITE(2,*)(ANORM(II),II=l,NN)       *PROVIDES   (ATPA)"1    IN  VECTOR  FORM 

CALL  MAT(ANORM.NMAX) 
C 

c  @@@@@@@@@@(a(^@@(a(a(a(a 

C  @  '  <a 

C  @   COMPUTE  VARIOUS  STATISTICAL  QUANTITIES  INCLUDING  THE     @ 

C  @    VELOCITY  OF  PEGASUS.  @ 

c   @  '  (a 

c 

c 

STOP 
END 


c 

c   @  @ 

C   @      SUBROUTINES  USED  IN  THE  PROGRAM      @ 
C   @  <a 

c 

SUBROUTINE  APPRO(TRANS ,X, Y , Z , SV.NT, AUX ,XP, YP, ZP) 

IMPLICIT  REAL*8(A-H,0-Z) 
C 

C  THIS  SUBROUTINE  COMPUTES  THE  APPROXIMATE  COORDINATES  FOR  PEGASUS 
C  INPUT:  TRANS  -  VECTOR  OF  TRANSPONDER  TIME  DELAYS  FOR  THE  DESIRED 
C  PEGASUS  POSITION. 
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C  X.Y.Z  -  ARRAYS  CONTAINING  THE  COORDINATES  OF  THE  TRANSPONDERS 

C  SV      VELOCITY  OF  SOUND  IN  WATER 

C  NT      NUMBER  OF  TRANSPONDERS 

C   OUTPUT:  XP.YP.ZP  -  COORDINATES  OF  PEGASUS 

C 

DIMENSION  TRANS(NT),  X(NT) , Y(NT) , Z(NT) , AUX( 2) 
C 

DO  10  I  =  l.NT 
TRANS(I)  =  TRANS(I)*SV 
10  CONTINUE 

DO  20  I  =  1,2 

AUX(I)  =  (TRANS(I+1)**2  -  ( ZP-Z( 1  +  1) )**2  -  X(I+1)**2  -  Y(I+1)**2) 
1  -  (TRANS(I)**2  -  (ZP-Z(I))**2  -  X(I)**2  -  Y(I)**2) 
20  CONTINUE 

YP  =  0.5D0*(AUX(1)/(X(1)-X(2))  -  AUX( 2)/(X( 2) -X( 3) ) )/(  (  Y(  1) -Y(  2)  ) 
1   /(X(l)-X(2))  -  (Y(2)-Y(3))/(X(2)-X(3))) 
XP  =  0.5D0*(AUX(1)/(Y(1)-Y(2))  -  AUX( 2)/( Y( 2) -Y( 3) ) )/(  (X(  1)  -X(  2)  ) 
1   /(Y(l)-Y(2))  -  (X(2)-X(3))/(Y(2)-Y(3))) 
DO  30  I  =  1, NT 
TRANS(I)  =  TRANS(I)/SV 
30  CONTINUE 
RETURN 
END 

SUBROUTI NE  SETUP(  X ,  Y ,  Z  ,  XP ,  YP  ,  ZP ,  S V ,  N ,  NT ,  NV ,  NT3  ,  NMAX ,  A ,  S  ,  ICOL ,  IROW , 
1  VALUE) 

IMPLICIT  REAL*8(A-H,0-Z) 
C 

C   THIS  MATRIX  SETS  UP  THE  A  AND  S  MATRICES  NEEDED  FOR  THE  BLOCK 
C    BY  BLOCK  FORMATION  OF  THE  NORMAL  MATRIX. 

C       INPUT:   X(I),  I  =  1,NT  X  COORDINATES  FOR  THE  TRANSPONDERS 
C  Y(I),  I  =  l.NT  Y      "       " 

C  Z(I),  I  =  l.NT  Z      "       " 

C  XP.YP.ZP,  APPROXIMATE  COORDINATES  FOR  PEGASUS. 

C  SV        VELOCITY  OF  SOUND  THROUGH  WATER. 

C  N       =  NTH  POSITION  OF  PEGASUS. 

C  NT,  NV,  NT3  SAME  AS  IN  THE  MAIN  PROGRAM. 

C  NMAX    =  ACTUAL  #  OF  PEGASUS  POSITIONS  IN  THIS  DROP 

C       OUTPUT:  A  =  THE  PORTION  OF  THE  A  MATRIX  CORRESPONDING  TO  THE 
C  NTH  PEGASUS  POSITION. 

C  S  =  AS  ABOVE,  BUT  THE  OUTER  BAND  OF  THE  A  MATRIX. 

C  VALUE  =  THE  NON  ZERO  //  IN  THE  A  MATRIX  CORRESPONDING 

C  TO  THE  COLUMN  //  HELD  IN  ICOL  AND  THE  ROW  //HELD 

C  IN  IROW. 

C 

DIMENSION  X(NT),  Y(NT) ,  Z(NT),  A(NT,3),  S(NT,NT3),  VALUE(NV) 

INTEGER*2  ICOL(NV),  IROW(NV) 
C 

NROW  =  (N-1)*NT 

NCOL  =  (N-l)*3 

NCOL1  =  NMAX*3 

ICOUNT  =  (N-1)*6*NT 

DO  20  I  =  l.NT  • 

J  =  1 

K  =  (I-l)*3  +  1 

DIST  =  DSQRT((XP  -  X(I))**2  +  (YP  -  Y(I))**2  +  (ZP  -  Z(I))**2) 
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A(I,J)  =  (XP  -  X(I))/(DIST*SV) 
ICOUNT  =  ICOUNT  +  1 
VALUE (I COUNT)  =  A(I,J) 
IROW( ICOUNT)  =  NROW  +  I 
ICOL( ICOUNT)  =  NCOL  +  J 
A(I,J+1)  =  (YP  -  Y(I))/(DIST*SV) 
ICOUNT  =  ICOUNT  +  1 
VALUE( ICOUNT)  =  A(I,J+1) 
IROW( ICOUNT)  =  NROW  +  I 
ICOL( ICOUNT)  =  NCOL  +  J+l 
A(I,J+2)  =  (ZP  -  Z(I))/(DIST*SV) 
ICOUNT  =  ICOUNT  +  1 
VALUE( ICOUNT)  =  A(I,J+2) 
IROW( ICOUNT)  =  NROW  +  I 
ICOL( ICOUNT)  =  NCOL  +  J+2 
S(I,K)  =  -A(I,J) 
ICOUNT  =  ICOUNT  +  1 
VALUE( ICOUNT)  =  S(I,K) 
IROW( ICOUNT)   =  NROW  +  I 
ICOL( ICOUNT)   =  NCOL1  +  K 
S(I,K+1)  =  -A(I,J+1) 
ICOUNT  =  ICOUNT  +  1 
VALUE( ICOUNT)  =  S(I,K+1) 
IROW( ICOUNT)  =  NROW  +  I 
ICOL( ICOUNT)  =  NCOL1  +  K+l 
S(I,K+2)  =  -A(I,J+2) 
ICOUNT  =  ICOUNT  +  1 
VALUE( ICOUNT)  =  S(I,K+2) 
IROW( ICOUNT)  =  NROW  +  I 
ICOL( ICOUNT)  =  NCOL1  +  K+2 
20  CONTINUE 

WRITE(2,100) 

100  FORMATCl'  ,///,15X,'A  MATRIX  AND  THEN  THE  S  MATRIX,  BLOCK  BY', 
1  '  BLOCK') 

WRITE(2,101)((A(I,J),J=1,3),I=1,4) 

101  F0RMAT('  ',  T2,4(2X,3(F10.  8,3X),/)) 
WRITE(2,102)((S(I,J),J=1,12),I=1,4) 

102  FORMATC  ',  T2,4(12(F10.  7,1X),/)) 
RETURN 

END 

SUBROUTINE  NORM(  A,S  ,P,L,N,NT,NT3  ,NN,NP,NMAX,ATP,ATPA,ATPS,ATPL, 
1  STP,STPS,STPL,ANORM,XHAT) 

IMPLICIT  REAL*8(A-H,0-Z) 
C 

C   THIS  SUBROUTINE  FORMS  THE  NORMAL  EQUATIONS,  BLOCK  BY  BLOCK 
C      INPUT:    A  MATRIX  FROM  SUBROUTINE  SETUP 

p  Q         II  II  II  II  I 

C  P  THE  WEIGHT  MATRIX  CORRESPONDING  TO  THIS  BLOCK 

C  L  THE  L  MATRIX 

C  N,  NT,  NT3,  AS  FOR  SUBROUTINE  SETUP 

C  NN,  NP,  NMAX,  AS  FOR  MAIN  PROGRAM 

C      OUTPUT:   A  NUMBER  OF  AUXILLARY  MATRICES  USED  IN  THE  MATRIX 

C  MULTIPLICATIONS,  IE.  ,  ATP, ATPA, ATPS , ATPL.STP.STPS , 

C  AND  STPL. 

C  MATRICES  SPS  AND  SPL  WHICH  COME  IN  VIA  THE  COMMON 
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C  BLOCK  ARE  UPDATED. 

C  ANORM  -  THE  NORMAL  MATRIX  (IN  ARRAY  FORM) 

C  XHAT  -  THE  ATPL  VECTOR 

C 

DIMENSION  A(NT,3),  S(NT,NT3) ,P(NT,NT) , ATP(NT, 3) ,  ATPA(3,3), 

1  ATPS(3,NT3),  ATPLC3.1),  STP(NT3,NT),  STPS(NT3 ,NT3) ,STPL(NT3  , 1)  , 

2  ANORM(NN),XHAT(NP) 
DOUBLE  PRECISION  L(NT,1) 
COMMON  SPS( 12, 12),  SPL(12,1) 

C 

C    PERFORM  THE  MATRIX  MULTIPLICATIONS 

C 

CALL  ATB(A,P,ATP,NT,3,NT) 

CALL  AB(ATP,A,ATPA,3,NT,3) 

CALL  AB(ATP,S,ATPS,3,NT,NT3) 

CALL  AB(ATP,L,ATPL,3,NT,1) 

CALL  ATB(S,P,STP,NT,NT3,NT) 

CALL  AB(STP,S,STPS,NT3,NT,NT3) 

CALL  AB(STP,L,STPL,NT3,NT,1) 

WRITE(2  100)((ATPA(I,J),J=1,3),ATPL(I,1),I=1,3) 
100  FORMAT( '  ',/,15X,  ' ATPA  BLOCK' , 10X, ' ATPL  VECTOR',//, 
1  3(2X,3(D10.  3,2X),  4X,D10.3,/)) 

WRITE(2.110)((STPS(I,J),J=1,12),I=1,12) 
110  FORMAT('  ',/,15X,  ' STPS  BLOCK*,//, 
1  T2,12(12(D10.3,1X),/)) 
C 

C    PLACE  ATPA  AND  ATPL  IN  THEIR  APPROPRIATE  POSITION  IN  EITHER  THE 
C   NORMAL  MATRIX  OR  THE  COLUMN  VECTOR 
C 

NCOL  =  (N-l)*3  +  1 

II  =  NCOL*(NCOL  +  l)/2 
ANORM(Il)  =  ATPA(l.l) 

III  =  II  +  NCOL 
ANORM(Ill)  =  ATPA(1,2) 
AN0RM(I11  +  1)  =  ATPA(2,2) 
1111  =  111  +  1  +  NCOL 
ANORM(Illl)  =  ATPA(1,3) 
ANORMCI111  +  1)  =  ATPA(2,3) 
AN0RM(I111  +  2)  =  ATPA(3,3) 

C 

DO  10  I  =  1,NT3 

INT  =  (NMAX*3)  +  I 

J  =  INT*(INT  -  l)/2  +  1  +  (N-l)*3 

ANORM(J)  =  ATPS(l.I) 

ANORM(J+l)  =  ATPS(2,I) 

ANORM(J+2)  =  ATPS(3,I) 
10  CONTINUE 
C 

DO  20  I  =  1,3 

XHAT(NCOL+I-l)  =ATPL(I,1) 
20  CONTINUE 
C 

C   UPDATE  SPS  AND  SPL 
C 

DO  30  I  =  1.NT3 
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SPL(I.l)  =  SPL(I.l)  +  STPL(I.l) 
DO  30  J  =  1.NT3 

SPS(I.J)  =  SPS(I.J)  +  STPS(I.J) 
30  CONTINUE 
RETURN 
END 

SUBROUTINE  ATB( A , D , R , L ,M ,N) 

IMPLICIT  REAL*8(A-H,0-Z) 
C 

C   THIS  SUBROUTINE  COMPUTES  THE  MATRIX  PRODUCT  A'B 
C    INPUT:    MATRIX  A  (L  X  M) 
C  MATRIX  B  (L  X  N) 

C   OUTPUT   MATRIX  R  (M  X  N) 
C 

DIMENSION  A(L,M),B(L,N),R(M,N) 

DO  10  I  =  1,M 

DO  10  J  =  1,N 

R(I,J)  =  0.  ODO 

DO  5  K  =  1  ,  L 

R(I,J)  =  R(I,J)  +  A(K,I)*B(K,J) 
5  CONTINUE 
10  CONTINUE 

RETURN 

END 

SUBROUTINE  AB( A , B ,R , L,M,N) 

IMPLICIT  REAL*8(A-H,0-Z) 
C 

C   FORM  THE  MATRIX  PRODUCT  R  =  AB. 
C   THE  MATRICES  A  AND  B  ARE  RETURNED  UNCHANGED 
C 

DIMENSION  A(L,M),B(M,N),R(L,N) 

DO  5  I  =  1,L 

DO  5  J  =  1,N 

R(I  ,J)  =  0.  ODO 

DO  5  K  =  1,M 

R(I,J)  =  R(I,J)  +  A(I,K)*B(K,J) 
5  CONTINUE 

RETURN 

END 

SUBROUTINE  RES ID(  VALUE  ,  ICOL , IROW , B IGP , LI ,  XHAT ,  NMAX  ,  NT ,  NV , NP  ,NOBS  , 
1  P,VTP,V,SIG,V1,SIG0) 
IMPLICIT  REAL*8(A-H,0-Z) 
C 

C    COMPUTES  THE  RESIDUALS  ON  ALL  THE  TRANSPONDER  TIME  DELAY 
C   OBSERVATIONS  AND  THE  A  POSTERIORI  VARIANCE  OF  UNIT  WEIGHT. 
C       INPUT:  VALUE,  ICOL, IROW,-  MATRICES  DERIVED  IN  SUBROUTINE  SETUP 
C  BIGP  -  THE  FULL  WEIGHT  MATRIX  FOR  ALL  OBSERVATIONS 

C  LI     THE  VECTOR  OF  "COMPUTED-OBSERVED"  OBSERVATIONS 

C  XHAT  -  THE  LEAST  SQUARES  SOLUTION  VECTOR 

C  NMAX, NT, NV.NP, NOBS  -  AS  IN  THE  MAIN  PROGRAM 

C        AUX1LLARY  MATRICES  :  P,  VTP ,  V.SIG 
C       OUTPUT:  VI    -  THE  VECTOR  OF  RESIDUALS 
C  SIGO  -  THE  A  POSTERIORI  VARIANCE  OF  UNIT  WEIGHT 
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DIMENSION  VALUE(NV) ,XHAT(NP) , Vl(NOBS) , BIGP(NT.NOBS) , P(NT.NT)  , 
1  VTP(1,NT),SIG(1,1),V(NT,1) 
DOUBLE  PRECISION  Ll(NOBS) 
INTEGER*2  ICOL(NV),  IROW(NV) 

DO  10  I  =  l.NOBS 
V1(I)  =  0.  ODO 
10  CONTINUE 
IC  =  1 

15  DO  20  I  =  l.NV 

IF(IROW(I).NE.  IC)  GO  TO  18 

16  Vl(IC)  =  Vl(IC)  +  VALUE(I)*XHAT(ICOL(I)) 
GO  TO  20 

18    Vl(IC)  =  Vl(IC)  +  Ll(IC) 
IC  =  IC  +  1 

IF(IC.  LE.NOBS)  GO  TO  16 
20  CONTINUE 

Vl(NOBS)  =  Vl(NOBS)  +  Ll(NOBS) 
C 

SIGO  =  0.  ODO 
IC  =  0 

DO  40  I  =  l.NMAX 
DO  30  J  =  1, NT 
V(J,1)  =  V1(IC*NT  +  J) 
DO  30  K  =  1,NT 
P(J,K)  =  BIGP(J,IC*NT  +  K) 
30    CONTINUE 

CALL  ATB(V,P,VTP,NT,1,NT) 
CALL  AB(VTP,V,SIG,1,NT,1) 
SIGO  =  SIGO  +  SIG(l.l) 
IC  =  IC  +  1 
40  CONTINUE 

SIGO  =  SIG0/FL0AT(N0BS-NMAX*3) 

RETURN 

END 

SUBROUTINE  P2L2( P2 , LB2 , XO , ANORM , L2 , XHAT , NP , NN ) 

IMPLICIT  REAL*8(A-H,0-Z) 
C 

C   THIS  SUBROUTINE  COMPUTES  THE  L2  MATRIX  AND  ADDS  IT  TO  THE  COLUMN 
C   VECTOR.  IN  ADDITION  IT  INCREMENTS  THE  DIAGONAL  ELEMENTS  OF  THE 
C   NORMAL  MATRIX  TO  INCLUDE  THE  INFLUENCE  OF  P2 
C 

DIMENSION  P2(NP),  XO(NP) , ANORM(NN) , XHAT( NP , 1) 

DOUBLE  PRECISION  L2(NP),  LB2(NP) 
C 

DO  10  I  =  l.NP 

L2(I)  =  XO(I)  -  LB2(I) 

INT  =  I*(I+l)/2 

ANORM(INT)    =  ANORM(INT)    +   P2(I) 

P2L  =   P2(I)*L2(I) 

XHAT(I.l)    =  XHAT(I.l)    +   P2L 
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CONTINUE 

RETURN 

END 


SUBROUTINE  MAT(C.NMAX) 
C   *NOT  PART  OF  J  HANNAH'S  ORIGINAL  PROGRAM.   ADDED  BY  M.  HASKELL. 
C   THIS  PROGRAM  READS  APPROPRIATE  VALUES  FROM  THE  VARIANCE-COVARIANCE 
C    MATRIX  FILE  (ANORM)  AND  PLACES  THEM  IN  THE  MATRICES  NEEDED  TO  COM- 
C    PUTE  VARIANCES  OF  THE  CURRENT  VELOCITIES. 

IMPLICIT  REAL*8  (A-H.O-Z) 

DOUBLE  PRECISION  C(  1) 
C      DIMENSION  C(444153) 

I(K)  =  K*(K+l)/2 


K 


K  =  1 

ROW  = 

WRITE 

WRITE 

WRITE 

DO  15 
WRITE 
WRITE 


(2,*)C(I(K)),  C(I(K+1)-1),  C(I(K+2)-2) 
(2,*)C(I(K+1)-1),  C(I(K+1)),  C(I(K+2)-l) 
(2,*)C(I(K+2)-2),  C(I(K+2)-l),  C(I(K+2)) 
J  =  l.NMAX+1 

(2,*)C(I(K+3)-3),  C(I(K+4)-4),  C(I(K+5)-5) 
(2,*)C(I(K+3)-2),  C(I(K+4)-3),  C(I(K+5)-4) 
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WRITE  (2,*)C(I(K+3)-l),  C(I(K+4)-2),  C(I(K+5)-3) 

WRITE  (2,*)C(I(K+3)),  C(I(K+4)-l),  C(I(K+5)-2) 
WRITE  (2,")C(I(K+4)-l) ,  C(I(K+4)),  C(l(K+5)-l) 
WRITE  (2,*)C(I(K+5)-2),  C(I(K+5)-l),  C(I(K+5)) 

K  =  K+3 

CONTINUE 


WRITE  (2,*)C(I(K+3)-6 
WRITE  (2,*)C(I(K+3)-5 
WRITE  (2,*)C(I(K+3)-4 

WRITE  (2,*)C(I(K+3)-3 
WRITE  (2,*)C(I(K+3)-2 
WRITE  (2,*)C(I(K+3)-l 

WRITE  (2,*)C(I(K+3)), 
WRITE  (2,*)C(I(K+4)-l 
WRITE  (2,*)C(I(K+5)-2 

K  =  K+3 

WRITE  (2,*)C(I(K+3)-9 
WRITE  (2,*)C(I(K+3)-8 
WRITE  (2,*)C(I(K+3)-7 

WRITE  (2,*)C(I(K+3)-6 
WRITE  (2,*)C(I(K+3)-5 
WRITE  (2,*)C(I(K+3)-4 

WRITE  (2,*)C(I(K+3)-3 
WRITE  (2,*)C(I(K+3)-2 
WRITE  (2,*)C(I(K+3)-l 


C(I(K+4)-l),  C(I(K+5)-2) 
C(I(K+4))(  C(I(K+5)-l) 
C(I(K+5)-l),  C(I(K+5)) 


C(I(K+4)-7),  C(I(K+5)-8) 
C(I(K+4)-6),  C(I(K+5)-7) 

C(I(K+4)-5),  C(I(K+5)-6) 


C(I(K+4) 

-A), 

C(I(K+5)-5) 

C(I(K+4) 

-3), 

C(I(K+5)-4) 

C(I(K+4) 

-2), 

C(I(K+5)-3) 

C(I(K+4)-10) 
C(I(K+4)-9), 

C(I(K+4)-8), 

C(I(K+4)-7), 
C(I(K+4)-6), 
C(I(K+4)-5), 

C(I(K+4)-4), 
C(I(K+4)-3), 
C(I(K+4)-2), 


,  C(I(K+5)-ll) 
C(I(K+5)-10) 
C(I(K+5)-9) 

C(I(K+5)-8) 
C(I(K+5)-7) 
C(I(K+5)-6) 

C(I(K+5)-5) 
C(I(K+5)-4) 
C(I(K+5)-3) 
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WRITE  (2,*)C(I(K+3)),  C(I(K+4)-1),  C(I(K+5)-2) 
WRITE  (2,*)C(I(K+4)-l),  C(I(K+4)),  C(I(K+5)-l) 
WRITE  (2,*)C(I(K+5)-2),  CCI(K+5)-l),  C(I(K+5)) 
C 

C      K  =  36 
C      DO  15  J  =  1,6 

C  WRITE  (2,*)C(I(K+3)-3),  C(I(K+4)-4),  C(I(K+5)-5) 
C  WRITE  (2,*)C(I(K+3)-2),  C(I(K+4)-3),  C(I(K+5)-A) 
C  WRITE  (2,*)C(I(K+3)-l),  C(I(K+4)-2),  C(I(K+5)-3) 
C 

C        WRITE  (2,*)C(I(K+3)),  C(I(K+4)-l),  C(I(K+5)-2) 
C        WRITE  (2,*)C(I(K+4)-l),  C(I(K+4)),  C(I(K+5)-l) 
C        WRITE  (2,*)C(I(K+5)-2),  C(I(K+5)-l),  C(I(K+5)) 
C      K  =  K+36 
C  15    CONTINUE 
RETURN 
END 
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APPENDIX  B 
CURRENT  VELOCITY  VARIANCE-COVARIANCE  PROGRAM 

This  program  computes  the  variance-covariance  matrix  of 
the  current  velocities  by  a  procedure  described  in  Chapter  VI, 
Section  D.  It  uses  data  from  the  Pegasus  position  variance- 
covariance  matrix  produced  by  the  Fortran  program  Pegasus  as 
described  in  Appendix  A. 
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C   THIS  PROGRAM  USES  DATA  FROM  THE  VARIANCE-COVARIANCE  MATRIX  PRODUCED 
C    BY  THE  PEGASUS  FORTRAN  PROGRAM.   IT  COMPUTES  THE  VARIANCE -COVAR I  - 
C    ANCE  MATRIX  OF  THE  CURRENT  VELOCITIES. 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL  SP1(3,3),  SP12(3,3),  SP2( 3 , 3) , V( 3, 3) , SU, SV.SW 
C 

N  =  1 

DO  100  I  =  1,3 

READ(1,*)(SP1(I,J),J  =  1,3) 
C         WRITE(2,120)(SP1(I,J),J  =  1,3) 

120     FORMAT(f  ' , 3X ,F10.  4 , 3X.F10.  4 , 3X.F10. 4) 
100  CONTINUE 
C 

10   DO  200  I  =  1,3 

READ(1,*,END=45)(SP12(I,J),J  =  1,3) 
C        WRITE(2,120)(SP12(I,J),J  =  1,3) 
200  CONTINUE 

DO  300  I  =  1,3 

READ(1,*)(SP2(I,J),J  =  1,3) 
C        WRITE(2,120)(SP2(I,J),J  =  1,3) 

300  CONTINUE 
C 
CC    WRITE(2,15)N 

15   FORMAT(   '/,5X,' VELOCITY  VARIANCE-COVARIANCE  MATRIX  FOR  ', 
1' POSITION   ,13,/) 
DO  400  I  =  1,3 

DO  350  J  =  1,3 
V(I,J)=((SP1(I,J)+SP2(I,J)-SP12(I,J)-SP12(J,I))/(16**2) 
350  CONTINUE 
C     WRITE(2,120)(V(I,J),J  =  1,3) 
400  CONTINUE 

SU  =  SQRT(V(1,1)) 

SV  =  SQRT(V(2,2)) 

SW  =  SQRT(V(3,3)) 

CC    WRITE(2  130)SU,  SV,  SW 

CC130  FORMAT( '  ',/, 3X, ' SIGMA  U  =' ,F10.  4 , 3X, ' SIGMA  V  =',F10.4,3X, 
CC    1  'SIGMA  W  =' ,F10.4/) 

WRITE(2,600)N,SU,SV,SW 
600  FORMAT( '  ' ,10X,'P0S  '  , 13 , 3X.F7.  4 , 3X.F7. 4 , 3X.F7. 4) 
DO  450  I  =  1,3 
,  DO  460  J  =  1,3 

SP1(I,J)  =  SP2(I,J) 
460   CONTINUE 
C       WRITE(2.140)(SP1(I,J),J=1,3) 
C  140   FORNATC'  ' , 3X.F10.  4 , 3X.F10.  4 , 3X.F10. 4) 
450  CONTINUE 
N  =  N  +1 
GO  TO  10 
45   CONTINUE 
STOP 
END 
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